From b0792a2fe27a51e32bb9f8f284b961116ad7eb76 Mon Sep 17 00:00:00 2001 From: Devroop Kar Date: Fri, 24 Jul 2026 12:47:20 -0400 Subject: [PATCH] Added black/flake8 check and formating changes --- .flake8 | 2 +- .github/workflows/python_format_check.yml | 31 ++ .../experiment_2_dataset_creation.py | 81 +-- .../experiment_3_dataset_creation.py | 63 ++- .../experiment_4_dataset_creation.py | 41 +- .../experiment_5_dataset_creation.py | 76 ++- .../experiment_8_dataset_creation.py | 361 +++++++++---- .../hamiltonian_generator.py | 75 ++- .../experiment_9_dataset_creation.py | 497 ++++++++++++------ .../hamiltonian_generator.py | 139 +++-- .../experiment_10_dataset_creation.py | 376 ++++++++----- .../hamiltonian_generator.py | 189 ++++--- .../experiment_11_dataset_creation.py | 109 ++-- .../experiment_12_dataset_creation.py | 84 +-- .../experiment_13_dataset_creation.py | 70 ++- website/build_site.py | 56 +- 16 files changed, 1506 insertions(+), 744 deletions(-) create mode 100644 .github/workflows/python_format_check.yml diff --git a/.flake8 b/.flake8 index bfccdeb..5c5faf8 100644 --- a/.flake8 +++ b/.flake8 @@ -1,4 +1,4 @@ [flake8] max-line-length = 120 ignore = E203,W503 -exclude = .bctds +exclude = .bctds,.venv,website diff --git a/.github/workflows/python_format_check.yml b/.github/workflows/python_format_check.yml new file mode 100644 index 0000000..9723bba --- /dev/null +++ b/.github/workflows/python_format_check.yml @@ -0,0 +1,31 @@ +# This workflow will install Python dependencies, run tests and lint with a variety of Python versions +# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions + +name: Flake8 Check + +on: + push: + branches: [ main ] + pull_request: + branches: [ main ] + +jobs: + python-format-check: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: [ "3.12" ] + + steps: + - uses: actions/checkout@v3.1.0 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v4.3.0 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -e . + - name: Lint with flake8 + run: | + flake8 . diff --git a/src/dataset_creation_scripts/02_pulse_width_sweep/experiment_2_dataset_creation.py b/src/dataset_creation_scripts/02_pulse_width_sweep/experiment_2_dataset_creation.py index 7c8f3e6..aef8cc8 100644 --- a/src/dataset_creation_scripts/02_pulse_width_sweep/experiment_2_dataset_creation.py +++ b/src/dataset_creation_scripts/02_pulse_width_sweep/experiment_2_dataset_creation.py @@ -160,7 +160,6 @@ import numpy as np from tqdm import tqdm - # ONE chip only, which is why there is no sample_id column. Recorded here so the # fact is not lost; it also goes into metadata.json. SAMPLE_NAME = "Si_111_native_oxide" @@ -178,7 +177,7 @@ AWG_CLOCK_MHZ = 9830.4 # pulse generator's clock -> 1 tick ~= 0.1017 ns ADC_CLOCK_MHZ = 552.96 # digitiser's clock -> 1 sample ~= 1.8084 ns -N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down +N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down N_PULSE_WIDTHS = len(PULSE_WIDTH_TICKS_LIST) # 42 N_FREQUENCIES = len(FREQUENCY_LIST_MHZ) # 500 @@ -477,19 +476,27 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): n_total_rows = N_PULSE_WIDTHS * N_FREQUENCIES * N_TIME_SAMPLES # 21,000,000 - print(f"\nFormat : long (one row per time sample)") + print("\nFormat : long (one row per time sample)") print(f"Sample : {SAMPLE_NAME} ({SAMPLE_MATERIAL})") - print(f"Pulse widths : {N_PULSE_WIDTHS} ({PULSE_WIDTHS_NS[0]:.3f}-" - f"{PULSE_WIDTHS_NS[-1]:.3f} ns, from {PULSE_WIDTH_TICKS_LIST[0]}-" - f"{PULSE_WIDTH_TICKS_LIST[-1]} ticks)") - print(f"Frequencies : {N_FREQUENCIES} ({FREQUENCY_LIST_MHZ[0]}-" - f"{FREQUENCY_LIST_MHZ[-1]} MHz, step 1)") + print( + f"Pulse widths : {N_PULSE_WIDTHS} ({PULSE_WIDTHS_NS[0]:.3f}-" + f"{PULSE_WIDTHS_NS[-1]:.3f} ns, from {PULSE_WIDTH_TICKS_LIST[0]}-" + f"{PULSE_WIDTH_TICKS_LIST[-1]} ticks)" + ) + print( + f"Frequencies : {N_FREQUENCIES} ({FREQUENCY_LIST_MHZ[0]}-" + f"{FREQUENCY_LIST_MHZ[-1]} MHz, step 1)" + ) print(f"Time samples : {N_TIME_SAMPLES} per measurement") - print(f"Total rows : {N_PULSE_WIDTHS} x {N_FREQUENCIES} x " - f"{N_TIME_SAMPLES} = {n_total_rows:,}") + print( + f"Total rows : {N_PULSE_WIDTHS} x {N_FREQUENCIES} x " + f"{N_TIME_SAMPLES} = {n_total_rows:,}" + ) print(f"Columns : {COLUMN_NAMES}") - print(f"Memory needed : ~{n_total_rows*len(COLUMN_NAMES)*4/1024**3:.2f} GB " - f"for the array alone\n") + print( + f"Memory needed : ~{n_total_rows * len(COLUMN_NAMES) * 4 / 1024**3:.2f} GB " + f"for the array alone\n" + ) # Make the whole output table up front and fill it in, rather than growing it # row by row (which would be far slower). float32 = 4 bytes per number, so @@ -522,8 +529,10 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): # truncated acquisition, or one of the DERIVED products that sit in the # same directory as the raw sweeps. if IQ_matrix.shape != (2, N_FREQUENCIES, N_TIME_SAMPLES): - print(f" WARNING: {filename} has shape {IQ_matrix.shape}, expected " - f"(2, {N_FREQUENCIES}, {N_TIME_SAMPLES}) — skipping") + print( + f" WARNING: {filename} has shape {IQ_matrix.shape}, expected " + f"(2, {N_FREQUENCIES}, {N_TIME_SAMPLES}) — skipping" + ) continue I_matrix = IQ_matrix[0] # shape (500, 1000) — I for every frequency @@ -575,9 +584,9 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): with open(cols_path, "w", encoding="utf-8") as fh: json.dump(build_column_info(), fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)}) — rows x columns") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") @@ -585,7 +594,7 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): # These do not change the data. They re-derive facts we already expect to be # true and shout if any of them is not, so that a silently broken conversion # (a missing file, a wrong axis, a rounding trap) cannot pass unnoticed. - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] # Did every file arrive? 42 files x 500 frequencies x 1000 instants should @@ -596,15 +605,17 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): f"— likely caused by missing .npy files" ) else: - print(f" [PASS] Row count: {row_idx:,} = {N_PULSE_WIDTHS} pulse widths x " - f"{N_FREQUENCIES} frequencies x {N_TIME_SAMPLES} samples") + print( + f" [PASS] Row count: {row_idx:,} = {N_PULSE_WIDTHS} pulse widths x " + f"{N_FREQUENCIES} frequencies x {N_TIME_SAMPLES} samples" + ) n_nan = int(np.sum(np.isnan(data_final))) n_inf = int(np.sum(np.isinf(data_final))) if n_nan > 0 or n_inf > 0: errors.append(f" [FAIL] Data contains {n_nan} NaN and {n_inf} Inf values") else: - print(f" [PASS] No NaN or Inf values in dataset") + print(" [PASS] No NaN or Inf values in dataset") # Are the measured numbers the right size? What comes BACK from the chip is a # few thousand a.u.; the 30000 number is how hard we SHOUT, not what we hear. @@ -620,8 +631,10 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): f"I=[{I_min:.0f}, {I_max:.0f}], Q=[{Q_min:.0f}, {Q_max:.0f}]" ) else: - print(f" [PASS] I/Q readout range: I=[{I_min:.1f}, {I_max:.1f}], " - f"Q=[{Q_min:.1f}, {Q_max:.1f}]") + print( + f" [PASS] I/Q readout range: I=[{I_min:.1f}, {I_max:.1f}], " + f"Q=[{Q_min:.1f}, {Q_max:.1f}]" + ) for col_idx, name, expected in [ (0, "pulse widths", N_PULSE_WIDTHS), @@ -642,8 +655,10 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): if n_at == 0: errors.append(" [FAIL] frequency_MHz == 4254 matches no rows after upcast") else: - print(f" [PASS] frequency_MHz exact under float64 upcast " - f"({n_at:,} rows at 4254 MHz)") + print( + f" [PASS] frequency_MHz exact under float64 upcast " + f"({n_at:,} rows at 4254 MHz)" + ) ts_col = data_final[:, 2] ts_min, ts_max = float(ts_col.min()), float(ts_col.max()) @@ -663,26 +678,30 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): cols_ok = list(reloaded["columns"]) == COLUMN_NAMES del reloaded if not same: - errors.append(f" [FAIL] Pickle round-trip: reloaded data differs from memory") + errors.append(" [FAIL] Pickle round-trip: reloaded data differs from memory") elif not cols_ok: - errors.append(f" [FAIL] Pickle round-trip: columns differ from memory") + errors.append(" [FAIL] Pickle round-trip: columns differ from memory") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly") + print(" [PASS] Pickle round-trip: reloaded array matches exactly") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ─────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") - print(f" (expect roughly {n_total_rows*55/1024**3:.1f} GB and a few minutes)") + print( + f" (expect roughly {n_total_rows * 55 / 1024**3:.1f} GB and a few minutes)" + ) write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) # ───────────────────────────────────────────────────────────────────────────── diff --git a/src/dataset_creation_scripts/03_amplitude_control/experiment_3_dataset_creation.py b/src/dataset_creation_scripts/03_amplitude_control/experiment_3_dataset_creation.py index 5260e80..7c5c7e2 100644 --- a/src/dataset_creation_scripts/03_amplitude_control/experiment_3_dataset_creation.py +++ b/src/dataset_creation_scripts/03_amplitude_control/experiment_3_dataset_creation.py @@ -198,7 +198,7 @@ PULSE_WIDTH_TICKS = 308 AWG_CLOCK_MHZ = 9830.4 # pulse generator's clock -> 1 tick ~= 0.1017 ns ADC_CLOCK_MHZ = 552.96 # digitiser's clock -> 1 sample ~= 1.8084 ns -N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down +N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down PULSE_WIDTH_NS = (PULSE_WIDTH_TICKS / AWG_CLOCK_MHZ) * 1000 # ~= 31.33 ns @@ -348,7 +348,10 @@ def build_column_info() -> dict: return { "dataset": "experiment_3_amplitude_control", "format": "long (one row per time sample)", - "n_rows": len(SAMPLES) * len(SPACING_LIST_NS) * len(AMPLITUDE_LIST) * N_TIME_SAMPLES, + "n_rows": len(SAMPLES) + * len(SPACING_LIST_NS) + * len(AMPLITUDE_LIST) + * N_TIME_SAMPLES, "n_columns": len(COLUMN_NAMES), "schema_rule": ( "The filename holds everything constant within a file (sample, pulse " @@ -527,14 +530,18 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): n_amps = len(AMPLITUDE_LIST) n_total_rows = n_samples * n_spacings * n_amps * N_TIME_SAMPLES # 1,800,000 - print(f"\nFormat : long (one row per time sample)") + print("\nFormat : long (one row per time sample)") print(f"Samples : {n_samples} ({SAMPLE_NAMES})") print(f"Spacings : {n_spacings} ({SPACING_LIST_NS} ns)") - print(f"Amplitudes : {n_amps} ({AMPLITUDE_LIST[0]}-{AMPLITUDE_LIST[-1]}, " - f"step 300)") + print( + f"Amplitudes : {n_amps} ({AMPLITUDE_LIST[0]}-{AMPLITUDE_LIST[-1]}, " + f"step 300)" + ) print(f"Time samples : {N_TIME_SAMPLES} per measurement") - print(f"Total rows : {n_samples} x {n_spacings} x {n_amps} x " - f"{N_TIME_SAMPLES} = {n_total_rows:,}") + print( + f"Total rows : {n_samples} x {n_spacings} x {n_amps} x " + f"{N_TIME_SAMPLES} = {n_total_rows:,}" + ) print(f"Columns : {COLUMN_NAMES}\n") # Make the whole output table up front and fill it in, rather than growing it @@ -638,9 +645,9 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): with open(cols_path, "w", encoding="utf-8") as fh: json.dump(build_column_info(), fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)}) — rows x columns") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") @@ -648,7 +655,7 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): # These do not change the data. They re-derive facts we already expect to be # true and shout if any of them is not, so that a silently broken conversion # (a missing file, a wrong axis, a rounding trap) cannot pass unnoticed. - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] # Did every file arrive? 18 files x 100 amplitudes x 1000 instants should @@ -659,15 +666,17 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): f"— likely caused by missing .npy files" ) else: - print(f" [PASS] Row count: {row_idx:,} = {n_samples} samples x {n_spacings} " - f"spacings x {n_amps} amplitudes x {N_TIME_SAMPLES} samples") + print( + f" [PASS] Row count: {row_idx:,} = {n_samples} samples x {n_spacings} " + f"spacings x {n_amps} amplitudes x {N_TIME_SAMPLES} samples" + ) n_nan = int(np.sum(np.isnan(data_final))) n_inf = int(np.sum(np.isinf(data_final))) if n_nan > 0 or n_inf > 0: errors.append(f" [FAIL] Data contains {n_nan} NaN and {n_inf} Inf values") else: - print(f" [PASS] No NaN or Inf values in dataset") + print(" [PASS] No NaN or Inf values in dataset") # Are the measured numbers the right size? What comes BACK from the chip is a # few thousand a.u.; the 30000 in amplitude_arb is how hard we SHOUT, not @@ -683,8 +692,10 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): f"I=[{I_min:.0f}, {I_max:.0f}], Q=[{Q_min:.0f}, {Q_max:.0f}]" ) else: - print(f" [PASS] I/Q readout range: I=[{I_min:.1f}, {I_max:.1f}], " - f"Q=[{Q_min:.1f}, {Q_max:.1f}]") + print( + f" [PASS] I/Q readout range: I=[{I_min:.1f}, {I_max:.1f}], " + f"Q=[{Q_min:.1f}, {Q_max:.1f}]" + ) for col_idx, name, expected in [ (COL["sample_id"], "sample IDs", n_samples), @@ -707,10 +718,12 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): spacing_64 = data_final[:, COL["spacing_ns"]].astype(np.float64) n_at_10 = int((spacing_64 == 10).sum()) if n_at_10 == 0: - errors.append(f" [FAIL] spacing_ns == 10 matches no rows after float64 upcast") + errors.append(" [FAIL] spacing_ns == 10 matches no rows after float64 upcast") else: - print(f" [PASS] spacing_ns exact under float64 upcast " - f"({n_at_10:,} rows at 10 ns)") + print( + f" [PASS] spacing_ns exact under float64 upcast " + f"({n_at_10:,} rows at 10 ns)" + ) ts_col = data_final[:, COL["timestamp_ns"]] ts_min, ts_max = float(ts_col.min()), float(ts_col.max()) @@ -727,25 +740,27 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): # — i.e. that nothing was lost or corrupted in the writing. reloaded = load_pickle(pkl_path) if not np.array_equal(reloaded["data"], data_final): - errors.append(f" [FAIL] Pickle round-trip: reloaded data differs from memory") + errors.append(" [FAIL] Pickle round-trip: reloaded data differs from memory") elif list(reloaded["columns"]) != COLUMN_NAMES: - errors.append(f" [FAIL] Pickle round-trip: columns differ from memory") + errors.append(" [FAIL] Pickle round-trip: columns differ from memory") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly") + print(" [PASS] Pickle round-trip: reloaded array matches exactly") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ─────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) # ───────────────────────────────────────────────────────────────────────────── diff --git a/src/dataset_creation_scripts/04_phase_control/experiment_4_dataset_creation.py b/src/dataset_creation_scripts/04_phase_control/experiment_4_dataset_creation.py index f72468b..5e8500d 100644 --- a/src/dataset_creation_scripts/04_phase_control/experiment_4_dataset_creation.py +++ b/src/dataset_creation_scripts/04_phase_control/experiment_4_dataset_creation.py @@ -181,7 +181,7 @@ PULSE_WIDTH_TICKS = 308 AWG_CLOCK_MHZ = 9830.4 # pulse generator's clock -> 1 tick ~= 0.1017 ns ADC_CLOCK_MHZ = 552.96 # digitiser's clock -> 1 sample ~= 1.8084 ns -N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down +N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down PULSE_WIDTH_NS = (PULSE_WIDTH_TICKS / AWG_CLOCK_MHZ) * 1000 # ~= 31.33 ns @@ -386,7 +386,10 @@ def build_column_info() -> dict: return { "dataset": "experiment_4_phase_control", "format": "long (one row per time sample)", - "n_rows": len(SAMPLES) * len(SPACING_LIST_NS) * len(PHASE_LIST_DEG) * N_TIME_SAMPLES, + "n_rows": len(SAMPLES) + * len(SPACING_LIST_NS) + * len(PHASE_LIST_DEG) + * N_TIME_SAMPLES, "n_columns": len(COLUMN_NAMES), "schema_rule": ( "The filename holds everything constant within a file (sample, pulse " @@ -577,10 +580,12 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): n_phases = len(PHASE_LIST_DEG) n_total_rows = n_samples * n_spacings * n_phases * N_TIME_SAMPLES # 2,178,000 - print(f"\nFormat : long (one row per time sample)") + print("\nFormat : long (one row per time sample)") print(f"Samples : {n_samples} ({SAMPLE_NAMES})") print(f"Spacings : {n_spacings} ({SPACING_LIST_NS} ns)") - print(f"Phases : {n_phases} ({PHASE_LIST_DEG[0]}-{PHASE_LIST_DEG[-1]} deg, step 3)") + print( + f"Phases : {n_phases} ({PHASE_LIST_DEG[0]}-{PHASE_LIST_DEG[-1]} deg, step 3)" + ) print(f"Time samples : {N_TIME_SAMPLES} per measurement") print( f"Total rows : {n_samples} x {n_spacings} x {n_phases} x " @@ -691,9 +696,9 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): with open(cols_path, "w", encoding="utf-8") as fh: json.dump(build_column_info(), fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)}) — rows x columns") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") @@ -701,7 +706,7 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): # These do not change the data. They re-derive facts we already expect to be # true and shout if any of them is not, so that a silently broken conversion # (a missing file, a wrong axis, a rounding trap) cannot pass unnoticed. - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] # Did every file arrive? 18 files x 121 phases x 1000 instants should give @@ -722,7 +727,7 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): if n_nan > 0 or n_inf > 0: errors.append(f" [FAIL] Data contains {n_nan} NaN and {n_inf} Inf values") else: - print(f" [PASS] No NaN or Inf values in dataset") + print(" [PASS] No NaN or Inf values in dataset") # Are the measured numbers the right size? What comes BACK from the chip is # a few thousand a.u.; the 30000 number elsewhere in this file is how hard @@ -764,9 +769,11 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): spacing_64 = data_final[:, COL["spacing_ns"]].astype(np.float64) n_at_10 = int((spacing_64 == 10).sum()) if n_at_10 == 0: - errors.append(f" [FAIL] spacing_ns == 10 matches no rows after float64 upcast") + errors.append(" [FAIL] spacing_ns == 10 matches no rows after float64 upcast") else: - print(f" [PASS] spacing_ns exact under float64 upcast ({n_at_10:,} rows at 10 ns)") + print( + f" [PASS] spacing_ns exact under float64 upcast ({n_at_10:,} rows at 10 ns)" + ) ts_col = data_final[:, COL["timestamp_ns"]] ts_min, ts_max = float(ts_col.min()), float(ts_col.max()) @@ -783,25 +790,27 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): # memory — i.e. that nothing was lost or corrupted in the writing. reloaded = load_pickle(pkl_path) if not np.array_equal(reloaded["data"], data_final): - errors.append(f" [FAIL] Pickle round-trip: reloaded data differs from memory") + errors.append(" [FAIL] Pickle round-trip: reloaded data differs from memory") elif list(reloaded["columns"]) != COLUMN_NAMES: - errors.append(f" [FAIL] Pickle round-trip: columns differ from memory") + errors.append(" [FAIL] Pickle round-trip: columns differ from memory") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly") + print(" [PASS] Pickle round-trip: reloaded array matches exactly") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ─────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) # ───────────────────────────────────────────────────────────────────────────── diff --git a/src/dataset_creation_scripts/05_phase_control_long_spacing/experiment_5_dataset_creation.py b/src/dataset_creation_scripts/05_phase_control_long_spacing/experiment_5_dataset_creation.py index 41a63b0..f4c1c17 100644 --- a/src/dataset_creation_scripts/05_phase_control_long_spacing/experiment_5_dataset_creation.py +++ b/src/dataset_creation_scripts/05_phase_control_long_spacing/experiment_5_dataset_creation.py @@ -209,7 +209,7 @@ PULSE_WIDTH_TICKS = 308 AWG_CLOCK_MHZ = 9830.4 # pulse generator's clock -> 1 tick ~= 0.1017 ns ADC_CLOCK_MHZ = 552.96 # digitiser's clock -> 1 sample ~= 1.8084 ns -N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down +N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down PULSE_WIDTH_NS = (PULSE_WIDTH_TICKS / AWG_CLOCK_MHZ) * 1000 # ~= 31.33 ns # The clock along each recorded trace. The digitiser samples every ~1.8 ns, so @@ -543,15 +543,17 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): n_phases = len(EXPECTED_PHASE_DEG) n_total_rows = n_spacings * n_phases * N_TIME_SAMPLES # 1,089,000 - print(f"\nFormat : long (one row per time sample)") + print("\nFormat : long (one row per time sample)") print(f"Sample : {SAMPLE_NAMES[0]} (sapphire + Shipley 1813)") print(f"Frequency : {DRIVE_FREQUENCY_MHZ} MHz") print(f"Temperature : {TEMPERATURE_MK} mK") print(f"Spacings : {n_spacings} ({SPACING_LIST_NS} ns)") print(f"Phases : {n_phases} (0-360 deg, step 3)") print(f"Time samples : {N_TIME_SAMPLES} per measurement") - print(f"Total rows : {n_spacings} x {n_phases} x {N_TIME_SAMPLES} " - f"= {n_total_rows:,}") + print( + f"Total rows : {n_spacings} x {n_phases} x {N_TIME_SAMPLES} " + f"= {n_total_rows:,}" + ) print(f"Columns : {COLUMN_NAMES}\n") # ── Pass 1: find the run start, so elapsed_s is relative to it ──────────── @@ -573,7 +575,7 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): # ── Pass 2: build ──────────────────────────────────────────────────────── data_all = np.empty((n_total_rows, len(COLUMN_NAMES)), dtype=np.float32) row_idx = 0 - sample_id = SAMPLES[0][1] + # sample_id = SAMPLES[0][1] for spacing_ns in tqdm(SPACING_LIST_NS, desc="Spacings"): @@ -590,14 +592,18 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): time_stamps = z["time_stamp_list"] if IQ_matrix.shape != (2, n_phases, N_TIME_SAMPLES): - print(f" WARNING: {filename} has shape {IQ_matrix.shape}, expected " - f"(2, {n_phases}, {N_TIME_SAMPLES}) — skipping") + print( + f" WARNING: {filename} has shape {IQ_matrix.shape}, expected " + f"(2, {n_phases}, {N_TIME_SAMPLES}) — skipping" + ) continue # The phase axis is in the file; verify rather than assume. if not np.array_equal(phase_array, EXPECTED_PHASE_DEG): - print(f" WARNING: {filename} pulse_phase_array does not match " - f"np.arange(0, 361, 3) — using the file's own axis") + print( + f" WARNING: {filename} pulse_phase_array does not match " + f"np.arange(0, 361, 3) — using the file's own axis" + ) # Wall-clock time, i.e. when this reading was actually taken in the lab. # Stored as SECONDS SINCE THE RUN STARTED rather than a calendar date: @@ -644,21 +650,23 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): with open(cols_path, "w", encoding="utf-8") as fh: json.dump(build_column_info(), fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)})") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") # ── Sanity checks ──────────────────────────────────────────────────────── - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] if row_idx != n_total_rows: errors.append(f" [FAIL] Row count: got {row_idx:,}, expected {n_total_rows:,}") else: - print(f" [PASS] Row count: {row_idx:,} = {n_spacings} spacings x " - f"{n_phases} phases x {N_TIME_SAMPLES} samples") + print( + f" [PASS] Row count: {row_idx:,} = {n_spacings} spacings x " + f"{n_phases} phases x {N_TIME_SAMPLES} samples" + ) # I and Q must never be NaN. (Unlike experiment 11, nothing in this # experiment is legitimately missing.) @@ -667,12 +675,18 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): if n_nan: errors.append(f" [FAIL] I/Q contain {n_nan} NaN values") else: - print(f" [PASS] No NaN in I/Q") + print(" [PASS] No NaN in I/Q") - I_min, I_max = float(data_final[:, COL["I"]].min()), float(data_final[:, COL["I"]].max()) - Q_min, Q_max = float(data_final[:, COL["Q"]].min()), float(data_final[:, COL["Q"]].max()) - print(f" [INFO] I range: [{I_min:.1f}, {I_max:.1f}] " - f"Q range: [{Q_min:.1f}, {Q_max:.1f}]") + I_min, I_max = float(data_final[:, COL["I"]].min()), float( + data_final[:, COL["I"]].max() + ) + Q_min, Q_max = float(data_final[:, COL["Q"]].min()), float( + data_final[:, COL["Q"]].max() + ) + print( + f" [INFO] I range: [{I_min:.1f}, {I_max:.1f}] " + f"Q range: [{Q_min:.1f}, {Q_max:.1f}]" + ) for col_idx, name, expected in [ (COL["spacing_ns"], "spacings", n_spacings), @@ -689,34 +703,40 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): spacing_64 = data_final[:, COL["spacing_ns"]].astype(np.float64) n_at_600 = int((spacing_64 == 600).sum()) if n_at_600 == 0: - errors.append(f" [FAIL] spacing_ns == 600 matches no rows after upcast") + errors.append(" [FAIL] spacing_ns == 600 matches no rows after upcast") else: - print(f" [PASS] spacing_ns exact under float64 upcast " - f"({n_at_600:,} rows at 600 ns)") + print( + f" [PASS] spacing_ns exact under float64 upcast " + f"({n_at_600:,} rows at 600 ns)" + ) el = data_final[:, COL["elapsed_s"]] - print(f" [INFO] elapsed_s range: [{el.min():.0f}, {el.max():.0f}] s " - f"({(el.max()-el.min())/60:.1f} min run)") + print( + f" [INFO] elapsed_s range: [{el.min():.0f}, {el.max():.0f}] s " + f"({(el.max() - el.min()) / 60:.1f} min run)" + ) reloaded = load_pickle(pkl_path) if not np.array_equal(reloaded["data"], data_final): errors.append(" [FAIL] Pickle round-trip: reloaded data differs") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly") + print(" [PASS] Pickle round-trip: reloaded array matches exactly") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) if __name__ == "__main__": diff --git a/src/dataset_creation_scripts/08_amplitude_control/experiment_8_dataset_creation.py b/src/dataset_creation_scripts/08_amplitude_control/experiment_8_dataset_creation.py index 502735c..bd8f232 100644 --- a/src/dataset_creation_scripts/08_amplitude_control/experiment_8_dataset_creation.py +++ b/src/dataset_creation_scripts/08_amplitude_control/experiment_8_dataset_creation.py @@ -168,32 +168,32 @@ # ───────────────────────────────────────────────────────────────────────────── # Pulse timing (nanoseconds). -PULSE_RING = 200 # single pulse used in the ring-down map (panel 0) -PULSE1_NS = PULSE_RING # pulse-1 duration (panels 1, 2) -GAP_NS = 100 # idle gap AFTER pulse 1 ends (panel 2) -PULSE2_NS = 200 # pulse-2 duration (panel 2) +PULSE_RING = 200 # single pulse used in the ring-down map (panel 0) +PULSE1_NS = PULSE_RING # pulse-1 duration (panels 1, 2) +GAP_NS = 100 # idle gap AFTER pulse 1 ends (panel 2) +PULSE2_NS = 200 # pulse-2 duration (panel 2) # Drive-frequency sweep (panel 0). Treated numerically as GHz == ns^-1. F_MIN, F_MAX = 3.0, 5.0 -N_FREQ = 400 +N_FREQ = 400 # Time grid: 0 to 1600 ns, 1000 points. T_MAX_NS, N_T = 1600, 1000 # The N = 4 coupled-TLS ensemble. The seed fixes both the TLS frequencies and # the coupling matrix; the draw ORDER below fixes the ensemble deterministically. -N_TLS = 4 -SEED = 2072025 -J_MIN, J_MAX = -0.05, 0.05 # XX coupling range (J in [-0.05, 0.05] GHz) +N_TLS = 4 +SEED = 2072025 +J_MIN, J_MAX = -0.05, 0.05 # XX coupling range (J in [-0.05, 0.05] GHz) # Dissipation: collective decay GAMMA, pure dephasing GAMMA_PHI (off here). GAMMA, GAMMA_PHI = 0.002, 0.0 # Drive amplitudes (dimensionless). -AMP_BASE = 0.10 # fixed amplitude for the ring-down map (panel 0) -AMP1_FIXED = 0.10 # fixed pulse-1 amplitude for the two-pulse sweep (panel 2) -N_AMP = 60 -AMP_LO, AMP_HI = 0.0, 0.10 # A1 sweep (panel 1) and A2 sweep (panel 2) +AMP_BASE = 0.10 # fixed amplitude for the ring-down map (panel 0) +AMP1_FIXED = 0.10 # fixed pulse-1 amplitude for the two-pulse sweep (panel 2) +N_AMP = 60 +AMP_LO, AMP_HI = 0.0, 0.10 # A1 sweep (panel 1) and A2 sweep (panel 2) # Panel bookkeeping. PANEL_NAMES = ["ring_down_map", "amp_sweep", "two_pulse_sweep"] @@ -222,6 +222,7 @@ # ENSEMBLE SETUP — the exact RNG stream that fixes the ensemble # ───────────────────────────────────────────────────────────────────────────── + def build_ensemble(smoke_test=False): """Return (init_freqs, H_int). @@ -245,6 +246,7 @@ def build_ensemble(smoke_test=False): return init_freqs, J from hamiltonian_generator import build_spin_spin_interactions_random_distribution + H_int = build_spin_spin_interactions_random_distribution( N_TLS, J_MIN, J_MAX, alpha_x=1.0, alpha_y=0.0, alpha_z=0.0 ) @@ -257,7 +259,7 @@ def coupling_matrix_from_seed(): rng_state = np.random.get_state() try: np.random.seed(SEED) - _ = np.random.uniform(F_MIN, F_MAX, N_TLS) # consume init_freqs draw + _ = np.random.uniform(F_MIN, F_MAX, N_TLS) # consume init_freqs draw J = np.random.uniform(J_MIN, J_MAX, size=(N_TLS, N_TLS)) J = np.tril(J, -1) + np.tril(J, -1).T finally: @@ -271,29 +273,33 @@ def coupling_matrix_from_seed(): # a machine without qutip. The real run never touches these functions. # ───────────────────────────────────────────────────────────────────────────── + def _synthetic_single(freq, amp, tlist, init_freqs, pulse_ns): """Cheap analytic surrogate for run_simulation_single_pulse (smoke test).""" detune = np.min(np.abs(np.asarray(init_freqs) - freq)) - lorentz = 1.0 / (1.0 + (detune / 0.15) ** 2) # near-resonance envelope - drive = amp * lorentz * np.sin( - np.pi * np.clip(tlist / max(pulse_ns, 1e-9), 0, 1) - ) ** 2 - tail = np.where(tlist > pulse_ns, - amp * lorentz * np.exp(-GAMMA * (tlist - pulse_ns)), 0.0) + lorentz = 1.0 / (1.0 + (detune / 0.15) ** 2) # near-resonance envelope + drive = ( + amp * lorentz * np.sin(np.pi * np.clip(tlist / max(pulse_ns, 1e-9), 0, 1)) ** 2 + ) + tail = np.where( + tlist > pulse_ns, amp * lorentz * np.exp(-GAMMA * (tlist - pulse_ns)), 0.0 + ) pop = np.where(tlist <= pulse_ns, drive, tail) return np.clip(pop, 0.0, None).astype(np.float64) -def _synthetic_double(freq, amp1, amp2, pulse1_ns, gap_ns, pulse2_ns, - tlist, init_freqs): +def _synthetic_double( + freq, amp1, amp2, pulse1_ns, gap_ns, pulse2_ns, tlist, init_freqs +): """Cheap analytic surrogate for run_simulation_double_pulse (smoke test).""" p1 = _synthetic_single(freq, amp1, tlist, init_freqs, pulse1_ns) t2_start = pulse1_ns + gap_ns shifted = tlist - t2_start - p2 = np.where(shifted >= 0, - _synthetic_single(freq, amp2, np.clip(shifted, 0, None), - init_freqs, pulse2_ns), - 0.0) + p2 = np.where( + shifted >= 0, + _synthetic_single(freq, amp2, np.clip(shifted, 0, None), init_freqs, pulse2_ns), + 0.0, + ) return np.clip(p1 * np.exp(-GAMMA * tlist) + p2, 0.0, None).astype(np.float64) @@ -301,8 +307,10 @@ def _synthetic_double(freq, amp1, amp2, pulse1_ns, gap_ns, pulse2_ns, # SWEEP RUNNERS — thin wrappers around the ORIGINAL physics functions # ───────────────────────────────────────────────────────────────────────────── -def _run_sweep(kind, param_list, freq_star, tlist, init_freqs, H_int, - workers, smoke_test): + +def _run_sweep( + kind, param_list, freq_star, tlist, init_freqs, H_int, workers, smoke_test +): """Run one panel's sweep, returning a (len(param_list), len(tlist)) array. kind = "freq" -> single-pulse ring-down map, param_list = drive freqs @@ -323,27 +331,64 @@ def _run_sweep(kind, param_list, freq_star, tlist, init_freqs, H_int, elif kind == "amp1": out[i] = _synthetic_single(freq_star, p, tlist, init_freqs, PULSE1_NS) else: # amp2 - out[i] = _synthetic_double(freq_star, AMP1_FIXED, p, - PULSE1_NS, GAP_NS, PULSE2_NS, - tlist, init_freqs) + out[i] = _synthetic_double( + freq_star, + AMP1_FIXED, + p, + PULSE1_NS, + GAP_NS, + PULSE2_NS, + tlist, + init_freqs, + ) return out # ── real physics: parallel qutip solver ────────────────────────────────── from concurrent.futures import ProcessPoolExecutor, as_completed from hamiltonian_generator import ( - run_simulation_single_pulse, run_simulation_double_pulse, + run_simulation_single_pulse, + run_simulation_double_pulse, ) def submit(pool, p): if kind == "freq": - return pool.submit(run_simulation_single_pulse, p, AMP_BASE, - tlist, init_freqs, H_int, GAMMA, GAMMA_PHI, PULSE1_NS) + return pool.submit( + run_simulation_single_pulse, + p, + AMP_BASE, + tlist, + init_freqs, + H_int, + GAMMA, + GAMMA_PHI, + PULSE1_NS, + ) if kind == "amp1": - return pool.submit(run_simulation_single_pulse, freq_star, p, - tlist, init_freqs, H_int, GAMMA, GAMMA_PHI, PULSE1_NS) - return pool.submit(run_simulation_double_pulse, freq_star, AMP1_FIXED, p, - PULSE1_NS, GAP_NS, PULSE2_NS, - tlist, init_freqs, H_int, GAMMA, GAMMA_PHI) + return pool.submit( + run_simulation_single_pulse, + freq_star, + p, + tlist, + init_freqs, + H_int, + GAMMA, + GAMMA_PHI, + PULSE1_NS, + ) + return pool.submit( + run_simulation_double_pulse, + freq_star, + AMP1_FIXED, + p, + PULSE1_NS, + GAP_NS, + PULSE2_NS, + tlist, + init_freqs, + H_int, + GAMMA, + GAMMA_PHI, + ) n_workers = max(1, min(workers, len(param_list))) with ProcessPoolExecutor(n_workers) as pool: @@ -357,8 +402,8 @@ def submit(pool, p): # METADATA # ───────────────────────────────────────────────────────────────────────────── -def build_metadata(freq_star, init_freqs, n_freq, n_amp, n_t, - smoke_test=False) -> dict: + +def build_metadata(freq_star, init_freqs, n_freq, n_amp, n_t, smoke_test=False) -> dict: """Instance metadata: the schema header, the resolved runtime values (actual TLS frequencies, FREQ_STAR, coupling matrix), and provenance. The full data-feature schema (grouped rows + colour legend) lives in @@ -420,7 +465,7 @@ def build_metadata(freq_star, init_freqs, n_freq, n_amp, n_t, "provenance": { "physics_module": "hamiltonian_generator.py (QuTiP-5 API shim only)", "attribution": "Uses code and simulation methods from the " - "Fitzpatrick Lab, Dartmouth College.", + "Fitzpatrick Lab, Dartmouth College.", "panel_c_note": "Panel c sweeps A2 (second pulse) with A1 fixed.", "no_new_physics": True, "smoke_test": smoke_test, @@ -461,71 +506,138 @@ def row(knob, dtype, nulls, code_var, value_range, in_file, notes, colour): { "section": "OVERVIEW", "rows": [ - row("what_it_shows", "Text", "-", "experiment_8_dataset_creation.py", - "Amplitude and multi-pulse control of a 4-TLS ensemble", "No", + row( + "what_it_shows", + "Text", + "-", + "experiment_8_dataset_creation.py", + "Amplitude and multi-pulse control of a 4-TLS ensemble", + "No", "a: ring-down map (single pulse); b: amplitude sweep at the " "optimal drive frequency (collapse-and-revival); c: two-pulse " - "protocol (A2 swept, A1 fixed, gap 100 ns).", "green"), + "protocol (A2 swept, A1 fixed, gap 100 ns).", + "green", + ), ], }, { "section": "PHYSICAL SYSTEM METADATA - the simulated TLS ensemble", "rows": [ - row("system", "String", "No", "N_TLS", - "N = 4 coupled TLS, Lindblad master equation", "No", - "Ensemble size.", "green"), - row("TLS_frequencies", "Float", "No", "init_freqs", + row( + "system", + "String", + "No", + "N_TLS", + "N = 4 coupled TLS, Lindblad master equation", + "No", + "Ensemble size.", + "green", + ), + row( + "TLS_frequencies", + "Float", + "No", + "init_freqs", "uniform[3.0, 5.0] GHz, seed 2072025 -> " - "{3.49, 3.17, 4.10, 4.19} GHz", "No", + "{3.49, 3.17, 4.10, 4.19} GHz", + "No", "Bare frequencies drawn with a fixed seed; the resulting " - "w_i/2pi = 3.49, 3.17, 4.10, 4.19 GHz.", "green"), - row("coupling_J", "Float", "No", + "w_i/2pi = 3.49, 3.17, 4.10, 4.19 GHz.", + "green", + ), + row( + "coupling_J", + "Float", + "No", "build_spin_spin_interactions_random_distribution", "J/2pi in [-50, 50] MHz (code uniform(-0.05, 0.05)), XX only", - "No", "Random symmetric dipole-dipole couplings between every " - "pair.", "green"), - row("dissipation_gamma", "Float", "No", "GAMMA, GAMMA_PHI", - "Gamma/2pi = 2.0 MHz (code 0.002); gamma_phi = 0", "No", - "Collective relaxation rate.", "green"), + "No", + "Random symmetric dipole-dipole couplings between every " + "pair.", + "green", + ), + row( + "dissipation_gamma", + "Float", + "No", + "GAMMA, GAMMA_PHI", + "Gamma/2pi = 2.0 MHz (code 0.002); gamma_phi = 0", + "No", + "Collective relaxation rate.", + "green", + ), ], }, { "section": "DRIVE & PULSE PARAMETERS", "rows": [ - row("drive_frequency", "Float", "No", "FREQ_AXIS / FREQ_STAR", + row( + "drive_frequency", + "Float", + "No", + "FREQ_AXIS / FREQ_STAR", "3.0-5.0 GHz, 400 pts (panel a); optimal wd/2pi = 4.15 GHz " - "(panels b,c)", "No", + "(panels b,c)", + "No", "Swept in the ring-down map; fixed at the ring-down-" - "maximizing freq for b,c.", "green"), - row("drive_amplitude", "Float", "No", + "maximizing freq for b,c.", + "green", + ), + row( + "drive_amplitude", + "Float", + "No", "AMP_BASE / AMP_SWEEP / AMP2_SWEEP", "A1/2pi = 100 MHz base (code 0.10); sweeps 0 -> 0.10 in 60 " - "steps", "No", + "steps", + "No", "Fixed for the ring-down map; swept 0->0.10 for the amplitude " - "panels (pulse-1 in b, pulse-2 in c).", "green"), - row("pulse_timing", "Float", "No", + "panels (pulse-1 in b, pulse-2 in c).", + "green", + ), + row( + "pulse_timing", + "Float", + "No", "PULSE_RING / PULSE1_NS / GAP_NS / PULSE2_NS", "single pulse: 200 ns; two-pulse: 200 ns / gap 100 ns / " - "200 ns", "No", + "200 ns", + "No", "tau = 200 ns; two-pulse protocol with inter-pulse gap " - "tau_g = 100 ns.", "green"), + "tau_g = 100 ns.", + "green", + ), ], }, { "section": "TIME GRID", "rows": [ - row("time_grid", "Float", "No", "T_MAX_NS, N_T", - "0 to 1600 ns, 1000 points", "No", - "Simulation time axis.", "blue"), + row( + "time_grid", + "Float", + "No", + "T_MAX_NS, N_T", + "0 to 1600 ns, 1000 points", + "No", + "Simulation time axis.", + "blue", + ), ], }, { "section": "COMPUTED OUTPUT ARRAY", "rows": [ - row("population", "Float", "No", "pop maps", + row( + "population", + "Float", + "No", + "pop maps", " vs (freq,time), (amp,time), (amp2,time)", - "No", "Collective excitation; ring-down, amplitude sweep, " - "two-pulse sweep.", "orange"), + "No", + "Collective excitation; ring-down, amplitude sweep, " + "two-pulse sweep.", + "orange", + ), ], }, ], @@ -559,6 +671,7 @@ def row(knob, dtype, nulls, code_var, value_range, in_file, notes, colour): # I/O # ───────────────────────────────────────────────────────────────────────────── + class _Tee: """Duplicate everything written to stdout into a log file as well, so each run leaves a persistent, human-readable record for reverification (the run @@ -596,7 +709,7 @@ def write_csv(data: np.ndarray, columns: list, path: str, chunk_rows: int = 500_ first = True with open(path, "w", newline="", encoding="utf-8") as fh: for start in tqdm(range(0, n, chunk_rows), desc="CSV", unit="chunk"): - df = pd.DataFrame(data[start: start + chunk_rows], columns=columns) + df = pd.DataFrame(data[start : start + chunk_rows], columns=columns) df.to_csv(fh, index=False, header=first) first = False @@ -605,6 +718,7 @@ def write_csv(data: np.ndarray, columns: list, path: str, chunk_rows: int = 500_ # MAIN # ───────────────────────────────────────────────────────────────────────────── + def build_dataset(output_dir, workers=4, write_csv_file=True, smoke_test=False): """Public entry point. Tees all console output to a run log for reverification, then delegates to _build_dataset_impl. The log is written to @@ -618,12 +732,15 @@ def build_dataset(output_dir, workers=4, write_csv_file=True, smoke_test=False): with open(log_path, "w", encoding="utf-8") as log_fh: sys.stdout = _Tee(original_stdout, log_fh) try: - print(f"# Experiment 8 run log — " - f"{datetime.datetime.now():%Y-%m-%d %H:%M:%S} " - f"({'SMOKE TEST' if smoke_test else 'full run'}, " - f"workers={workers})") - result = _build_dataset_impl(output_dir, workers, write_csv_file, - smoke_test) + print( + f"# Experiment 8 run log — " + f"{datetime.datetime.now():%Y-%m-%d %H:%M:%S} " + f"({'SMOKE TEST' if smoke_test else 'full run'}, " + f"workers={workers})" + ) + result = _build_dataset_impl( + output_dir, workers, write_csv_file, smoke_test + ) finally: sys.stdout = original_stdout @@ -647,14 +764,14 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F freq_axis = np.linspace(F_MIN, F_MAX, n_freq) amp_sweep = np.linspace(AMP_LO, AMP_HI, n_amp) - tlist = np.linspace(0.0, T_MAX_NS, n_t) + tlist = np.linspace(0.0, T_MAX_NS, n_t) n_total_rows = (n_freq + n_amp + n_amp) * n_t print(f"\nExperiment 8 dataset builder {'(SMOKE TEST)' if smoke_test else ''}") - print(f"Format : long (one row per time sample)") + print("Format : long (one row per time sample)") print(f"System : N={N_TLS} coupled TLS, Lindblad master equation") - print(f"Panels : 0=ring-down map, 1=amp sweep, 2=two-pulse sweep") + print("Panels : 0=ring-down map, 1=amp sweep, 2=two-pulse sweep") print(f"Drive freqs : {n_freq} ({F_MIN}-{F_MAX} GHz, panel 0)") print(f"Amplitudes : {n_amp} ({AMP_LO}-{AMP_HI}, panels 1 & 2)") print(f"Time samples : {n_t} (0-{T_MAX_NS} ns)") @@ -667,24 +784,27 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F # ── panel 0: ring-down map + FREQ_STAR ─────────────────────────────────── print("[panel 0] ring-down map ...") - pop_ring = _run_sweep("freq", freq_axis, None, tlist, init_freqs, H_int, - workers, smoke_test) + pop_ring = _run_sweep( + "freq", freq_axis, None, tlist, init_freqs, H_int, workers, smoke_test + ) mask_tail = tlist > PULSE1_NS if mask_tail.any(): freq_star = float(freq_axis[np.argmax(pop_ring[:, mask_tail].sum(axis=1))]) - else: # tiny smoke grid may have no tail sample + else: # tiny smoke grid may have no tail sample freq_star = float(freq_axis[np.argmax(pop_ring.sum(axis=1))]) print(f" FREQ_STAR : {freq_star:.4f} GHz") # ── panel 1: single-pulse amplitude sweep ──────────────────────────────── print("\n[panel 1] amplitude sweep ...") - pop_amp = _run_sweep("amp1", amp_sweep, freq_star, tlist, init_freqs, H_int, - workers, smoke_test) + pop_amp = _run_sweep( + "amp1", amp_sweep, freq_star, tlist, init_freqs, H_int, workers, smoke_test + ) # ── panel 2: two-pulse (A2) sweep ──────────────────────────────────────── print("\n[panel 2] two-pulse sweep ...") - pop_double = _run_sweep("amp2", amp_sweep, freq_star, tlist, init_freqs, H_int, - workers, smoke_test) + pop_double = _run_sweep( + "amp2", amp_sweep, freq_star, tlist, init_freqs, H_int, workers, smoke_test + ) # ── assemble the long table ────────────────────────────────────────────── data_all = np.empty((n_total_rows, len(COLUMN_NAMES)), dtype=np.float32) @@ -696,7 +816,8 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, pop_matrix, param_iter): end = row_idx + n_t data_all[row_idx:end, COL["panel_id"]] = panel_id data_all[row_idx:end, COL["drive_frequency_GHz"]] = ( - drive_freq[row_i] if np.ndim(drive_freq) else drive_freq) + drive_freq[row_i] if np.ndim(drive_freq) else drive_freq + ) data_all[row_idx:end, COL["amp1"]] = a1[row_i] if np.ndim(a1) else a1 data_all[row_idx:end, COL["amp2"]] = a2[row_i] if np.ndim(a2) else a2 data_all[row_idx:end, COL["pulse1_ns"]] = p1 @@ -711,8 +832,17 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, pop_matrix, param_iter): # panel 1: amp1 swept at freq_star; single pulse emit(1, freq_star, amp_sweep, 0.0, PULSE1_NS, 0.0, 0.0, pop_amp, amp_sweep) # panel 2: amp2 swept at freq_star; two pulses, amp1 fixed - emit(2, freq_star, AMP1_FIXED, amp_sweep, PULSE1_NS, GAP_NS, PULSE2_NS, - pop_double, amp_sweep) + emit( + 2, + freq_star, + AMP1_FIXED, + amp_sweep, + PULSE1_NS, + GAP_NS, + PULSE2_NS, + pop_double, + amp_sweep, + ) data_final = data_all[:row_idx] @@ -736,14 +866,14 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, pop_matrix, param_iter): with open(cols_path, "w", encoding="utf-8") as fh: json.dump(column_doc, fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)})") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") # ── sanity checks ──────────────────────────────────────────────────────── - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] if row_idx != n_total_rows: @@ -756,7 +886,7 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, pop_matrix, param_iter): if n_nan or n_inf: errors.append(f" [FAIL] {n_nan} NaN and {n_inf} Inf values present") else: - print(f" [PASS] No NaN or Inf values") + print(" [PASS] No NaN or Inf values") pop = data_final[:, COL["population"]] if pop.min() < -1e-3: @@ -772,17 +902,21 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, pop_matrix, param_iter): ts = data_final[:, COL["timestamp_ns"]] if ts.min() < -1e-6 or abs(ts.max() - tlist[-1]) > 1e-3: - errors.append(f" [FAIL] Timestamp range: [{ts.min():.3f}, {ts.max():.3f}] " - f"(expected [0, {tlist[-1]:.3f}])") + errors.append( + f" [FAIL] Timestamp range: [{ts.min():.3f}, {ts.max():.3f}] " + f"(expected [0, {tlist[-1]:.3f}])" + ) else: print(f" [PASS] Timestamp range: [{ts.min():.3f}, {ts.max():.3f}] ns") # pulse2_ns and gap_ns must be non-zero only for panel 2. c_mask = data_final[:, COL["panel_id"]] == 2 - ok_p2 = (np.all(data_final[c_mask, COL["pulse2_ns"]] == PULSE2_NS) and - np.all(data_final[~c_mask, COL["pulse2_ns"]] == 0.0)) - ok_gap = (np.all(data_final[c_mask, COL["gap_ns"]] == GAP_NS) and - np.all(data_final[~c_mask, COL["gap_ns"]] == 0.0)) + ok_p2 = np.all(data_final[c_mask, COL["pulse2_ns"]] == PULSE2_NS) and np.all( + data_final[~c_mask, COL["pulse2_ns"]] == 0.0 + ) + ok_gap = np.all(data_final[c_mask, COL["gap_ns"]] == GAP_NS) and np.all( + data_final[~c_mask, COL["gap_ns"]] == 0.0 + ) if not (ok_p2 and ok_gap): errors.append(" [FAIL] gap_ns/pulse2_ns should be non-zero only for panel 2") else: @@ -798,21 +932,23 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, pop_matrix, param_iter): if not np.array_equal(reloaded["data"], data_final): errors.append(" [FAIL] Pickle round-trip: reloaded data differs") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly") + print(" [PASS] Pickle round-trip: reloaded array matches exactly") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) return pkl_path @@ -826,18 +962,25 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, pop_matrix, param_iter): "pickle and a CSV (N=4 TLS amplitude and multi-pulse control)." ) parser.add_argument( - "--output_dir", type=str, default=_script_dir, + "--output_dir", + type=str, + default=_script_dir, help="Where to write experiment_8/ (default: next to this script).", ) parser.add_argument( - "--workers", type=int, default=4, + "--workers", + type=int, + default=4, help="Parallel solver processes (real run only; ignored in smoke test).", ) parser.add_argument( - "--no_csv", action="store_true", help="Skip the CSV, write only the pickle.", + "--no_csv", + action="store_true", + help="Skip the CSV, write only the pickle.", ) parser.add_argument( - "--smoke-test", action="store_true", + "--smoke-test", + action="store_true", help="Tiny grids + synthetic population; no qutip required.", ) args = parser.parse_args() diff --git a/src/dataset_creation_scripts/08_amplitude_control/hamiltonian_generator.py b/src/dataset_creation_scripts/08_amplitude_control/hamiltonian_generator.py index 3d9b0c6..59f9cb2 100644 --- a/src/dataset_creation_scripts/08_amplitude_control/hamiltonian_generator.py +++ b/src/dataset_creation_scripts/08_amplitude_control/hamiltonian_generator.py @@ -3,18 +3,6 @@ ──────────────────────── QuTiP helpers for driven N-TLS simulations (Python ≥ 3.9, QuTiP ≥ 5). -Changes in this revision -──────────────────────── -• `run_simulation_double_pulse` takes - pulse1_ns, gap_ns, pulse2_ns (instead of pulse1_ns, delay_ns, pulse2_ns) - where `gap_ns` = time waited *after pulse 1 finishes*. -• Drive functions use the same envelope as the single-pulse definition. -""" -""" -hamiltonian_generator.py -──────────────────────── -QuTiP helpers for driven N-TLS simulations (Python ≥ 3.9, QuTiP ≥ 5). - Changes in this revision ──────────────────────── • `run_simulation_double_pulse` takes @@ -46,8 +34,17 @@ def T(op, k): # ────────────────────────── single pulse ───────────────────────────── -def run_simulation_single_pulse(freq, amp, tlist, init_freqs, interactions, - gamma, gamma_phi, pulse_ns, T_ramp_ns=1.0): +def run_simulation_single_pulse( + freq, + amp, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + pulse_ns, + T_ramp_ns=1.0, +): """ freq : ns⁻¹ (treat GHz numerically as ns⁻¹, so no 2π factor here) amp : drive amplitude @@ -59,6 +56,7 @@ def run_simulation_single_pulse(freq, amp, tlist, init_freqs, interactions, H_tls = sum(init_freqs[j] * sz[j] / 2 for j in range(N)) T_ramp = T_ramp_ns + def drive_coeff(t, _=None): if t < T_ramp: env = 0.5 * (1 - np.cos(np.pi * t / T_ramp)) @@ -85,10 +83,20 @@ def drive_coeff(t, _=None): # ───────────────────────── two pulses ──────────────────────────────── -def run_simulation_double_pulse(freq, amp1, amp2, - pulse1_ns, gap_ns, pulse2_ns, - tlist, init_freqs, interactions, - gamma, gamma_phi, T_ramp_ns=1.0): +def run_simulation_double_pulse( + freq, + amp1, + amp2, + pulse1_ns, + gap_ns, + pulse2_ns, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + T_ramp_ns=1.0, +): """ Two pulses at the same frequency. • Pulse 1: starts at t=0, duration = pulse1_ns @@ -100,9 +108,10 @@ def run_simulation_double_pulse(freq, amp1, amp2, Sx, Sm = sum(sx), sum(sm) H_tls = sum(init_freqs[j] * sz[j] / 2 for j in range(N)) - t_start_p2 = pulse1_ns + gap_ns # absolute start of pulse 2 + t_start_p2 = pulse1_ns + gap_ns # absolute start of pulse 2 T_ramp = T_ramp_ns + def drive_coeff(t, _=None): # ------- pulse 1 ------- if t < T_ramp: @@ -142,8 +151,9 @@ def drive_coeff(t, _=None): # ─────────────────── random XX couplings (αy,z optional) ───────────── -def build_spin_spin_interactions_random_distribution(N, J_min, J_max, - alpha_x=1.0, alpha_y=0.0, alpha_z=0.0): +def build_spin_spin_interactions_random_distribution( + N, J_min, J_max, alpha_x=1.0, alpha_y=0.0, alpha_z=0.0 +): _, sx, sy, sz = _tls_ops(N) J = np.random.uniform(J_min, J_max, size=(N, N)) J = np.tril(J, -1) + np.tril(J, -1).T @@ -153,15 +163,24 @@ def build_spin_spin_interactions_random_distribution(N, J_min, J_max, for j in range(i + 1, N): Jij = J[i, j] H_int += Jij * ( - alpha_x * sx[i] * sx[j] + - alpha_y * sy[i] * sy[j] + - alpha_z * sz[i] * sz[j] + alpha_x * sx[i] * sx[j] + + alpha_y * sy[i] * sy[j] + + alpha_z * sz[i] * sz[j] ) return H_int # ───────── legacy wrapper (unchanged signature) ───────── -def run_simulation_for_frequency(freq, tlist, init_freqs, interactions, - gamma, gamma_phi, drive_ampl, pulse_duration): - return run_simulation_single_pulse(freq, drive_ampl, tlist, init_freqs, - interactions, gamma, gamma_phi, pulse_duration) +def run_simulation_for_frequency( + freq, tlist, init_freqs, interactions, gamma, gamma_phi, drive_ampl, pulse_duration +): + return run_simulation_single_pulse( + freq, + drive_ampl, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + pulse_duration, + ) diff --git a/src/dataset_creation_scripts/09_separation_phase_sweep/experiment_9_dataset_creation.py b/src/dataset_creation_scripts/09_separation_phase_sweep/experiment_9_dataset_creation.py index dc40138..1129c7c 100644 --- a/src/dataset_creation_scripts/09_separation_phase_sweep/experiment_9_dataset_creation.py +++ b/src/dataset_creation_scripts/09_separation_phase_sweep/experiment_9_dataset_creation.py @@ -173,18 +173,18 @@ # ───────────────────────────────────────────────────────────────────────────── # Pulse timing (nanoseconds). Both pulses share one duration in this experiment. -PULSE_NS = 200 -FIXED_GAP_NS = 150 # gap held fixed while the phase is swept (panel 2) +PULSE_NS = 200 +FIXED_GAP_NS = 150 # gap held fixed while the phase is swept (panel 2) # Drive-frequency sweep (panel 0). Treated numerically as GHz == ns^-1. F_MIN, F_MAX = 2.0, 5.0 -N_FREQ = 400 +N_FREQ = 400 # Gap sweep (panel 1) and phase sweep (panel 2). GAP_LO, GAP_HI = 0.0, 800.0 -N_GAP = 120 +N_GAP = 120 PHASE_LO, PHASE_HI = 0.0, 2.0 * np.pi -N_PHASE = 120 +N_PHASE = 120 # Time grid: 0 to 1600 ns, 1000 points. T_MAX_NS, N_T = 1600, 1000 @@ -193,14 +193,14 @@ # [3.0, 4.0] after a throwaway random draw; the coupling is then drawn from the # seeded RNG. The draw ORDER (one 2-value uniform consumed first) is fixed below # so the coupling is deterministic on every run. -N_TLS = 2 -SEED = 42 -INIT_FREQS = [3.0, 4.0] -J_MIN, J_MAX = -0.05, 0.05 # XX coupling range +N_TLS = 2 +SEED = 42 +INIT_FREQS = [3.0, 4.0] +J_MIN, J_MAX = -0.05, 0.05 # XX coupling range # Dissipation and drive amplitude. GAMMA, GAMMA_PHI = 0.002, 0.0 -AMP_DRIVE = 0.12 # both pulses use this amplitude +AMP_DRIVE = 0.12 # both pulses use this amplitude # Panel bookkeeping. PANEL_NAMES = ["ring_down_map", "gap_sweep", "phase_sweep"] @@ -228,6 +228,7 @@ # ENSEMBLE SETUP — the exact RNG stream that fixes the ensemble # ───────────────────────────────────────────────────────────────────────────── + def build_ensemble(smoke_test=False): """Return (init_freqs, H_int). @@ -241,7 +242,7 @@ def build_ensemble(smoke_test=False): symmetric coupling matrix (a plain ndarray) instead of a QuTiP operator. """ np.random.seed(SEED) - _ = np.random.uniform(F_MIN, F_MAX, N_TLS) # consume the throwaway draw + _ = np.random.uniform(F_MIN, F_MAX, N_TLS) # consume the throwaway draw init_freqs = list(INIT_FREQS) if smoke_test: @@ -250,6 +251,7 @@ def build_ensemble(smoke_test=False): return init_freqs, J from hamiltonian_generator import build_spin_spin_interactions_random_distribution + H_int = build_spin_spin_interactions_random_distribution(N_TLS, J_MIN, J_MAX) return init_freqs, H_int @@ -260,7 +262,7 @@ def coupling_matrix_from_seed(): rng_state = np.random.get_state() try: np.random.seed(SEED) - _ = np.random.uniform(F_MIN, F_MAX, N_TLS) # consume throwaway draw + _ = np.random.uniform(F_MIN, F_MAX, N_TLS) # consume throwaway draw J = np.random.uniform(J_MIN, J_MAX, size=(N_TLS, N_TLS)) J = np.tril(J, -1) + np.tril(J, -1).T finally: @@ -273,30 +275,37 @@ def coupling_matrix_from_seed(): # NOT physically meaningful; it only exercises the dataset plumbing without qutip. # ───────────────────────────────────────────────────────────────────────────── + def _synthetic_single(freq, amp, tlist, init_freqs, pulse_ns): detune = np.min(np.abs(np.asarray(init_freqs) - freq)) lorentz = 1.0 / (1.0 + (detune / 0.15) ** 2) - drive = amp * lorentz * np.sin( - np.pi * np.clip(tlist / max(pulse_ns, 1e-9), 0, 1) - ) ** 2 - tail = np.where(tlist > pulse_ns, - amp * lorentz * np.exp(-GAMMA * (tlist - pulse_ns)), 0.0) + drive = ( + amp * lorentz * np.sin(np.pi * np.clip(tlist / max(pulse_ns, 1e-9), 0, 1)) ** 2 + ) + tail = np.where( + tlist > pulse_ns, amp * lorentz * np.exp(-GAMMA * (tlist - pulse_ns)), 0.0 + ) pop = np.where(tlist <= pulse_ns, drive, tail) return np.clip(pop, 0.0, None).astype(np.float64) -def _synthetic_double_phase(freq, amp1, amp2, pulse1_ns, gap_ns, pulse2_ns, - phase2_rad, tlist, init_freqs): +def _synthetic_double_phase( + freq, amp1, amp2, pulse1_ns, gap_ns, pulse2_ns, phase2_rad, tlist, init_freqs +): """Phase-aware surrogate: the second pulse interferes with the first with a (1 + cos(phase2))/2 weighting, so the phase axis actually varies the output.""" p1 = _synthetic_single(freq, amp1, tlist, init_freqs, pulse1_ns) t2_start = pulse1_ns + gap_ns shifted = tlist - t2_start interf = 0.5 * (1.0 + np.cos(phase2_rad)) - p2 = np.where(shifted >= 0, - interf * _synthetic_single(freq, amp2, np.clip(shifted, 0, None), - init_freqs, pulse2_ns), - 0.0) + p2 = np.where( + shifted >= 0, + interf + * _synthetic_single( + freq, amp2, np.clip(shifted, 0, None), init_freqs, pulse2_ns + ), + 0.0, + ) return np.clip(p1 * np.exp(-GAMMA * tlist) + p2, 0.0, None).astype(np.float64) @@ -304,8 +313,10 @@ def _synthetic_double_phase(freq, amp1, amp2, pulse1_ns, gap_ns, pulse2_ns, # SWEEP RUNNERS — thin wrappers around the ORIGINAL physics functions # ───────────────────────────────────────────────────────────────────────────── -def _run_sweep(kind, param_list, freq_star, tlist, init_freqs, H_int, - workers, smoke_test): + +def _run_sweep( + kind, param_list, freq_star, tlist, init_freqs, H_int, workers, smoke_test +): """Run one panel's sweep, returning a (len(param_list), len(tlist)) array. kind = "freq" -> single-pulse ring-down map, param_list = drive freqs @@ -322,33 +333,81 @@ def _run_sweep(kind, param_list, freq_star, tlist, init_freqs, H_int, if kind == "freq": out[i] = _synthetic_single(p, AMP_DRIVE, tlist, init_freqs, PULSE_NS) elif kind == "gap": - out[i] = _synthetic_double_phase(freq_star, AMP_DRIVE, AMP_DRIVE, - PULSE_NS, p, PULSE_NS, 0.0, - tlist, init_freqs) + out[i] = _synthetic_double_phase( + freq_star, + AMP_DRIVE, + AMP_DRIVE, + PULSE_NS, + p, + PULSE_NS, + 0.0, + tlist, + init_freqs, + ) else: # phase - out[i] = _synthetic_double_phase(freq_star, AMP_DRIVE, AMP_DRIVE, - PULSE_NS, FIXED_GAP_NS, PULSE_NS, p, - tlist, init_freqs) + out[i] = _synthetic_double_phase( + freq_star, + AMP_DRIVE, + AMP_DRIVE, + PULSE_NS, + FIXED_GAP_NS, + PULSE_NS, + p, + tlist, + init_freqs, + ) return out from concurrent.futures import ProcessPoolExecutor, as_completed from hamiltonian_generator import ( - run_simulation_single_pulse, run_simulation_double_pulse_phase, + run_simulation_single_pulse, + run_simulation_double_pulse_phase, ) def submit(pool, p): if kind == "freq": - return pool.submit(run_simulation_single_pulse, p, AMP_DRIVE, - tlist, init_freqs, H_int, GAMMA, GAMMA_PHI, PULSE_NS) + return pool.submit( + run_simulation_single_pulse, + p, + AMP_DRIVE, + tlist, + init_freqs, + H_int, + GAMMA, + GAMMA_PHI, + PULSE_NS, + ) if kind == "gap": - return pool.submit(run_simulation_double_pulse_phase, - freq_star, AMP_DRIVE, AMP_DRIVE, - PULSE_NS, p, PULSE_NS, 0.0, - tlist, init_freqs, H_int, GAMMA, GAMMA_PHI) - return pool.submit(run_simulation_double_pulse_phase, - freq_star, AMP_DRIVE, AMP_DRIVE, - PULSE_NS, FIXED_GAP_NS, PULSE_NS, p, - tlist, init_freqs, H_int, GAMMA, GAMMA_PHI) + return pool.submit( + run_simulation_double_pulse_phase, + freq_star, + AMP_DRIVE, + AMP_DRIVE, + PULSE_NS, + p, + PULSE_NS, + 0.0, + tlist, + init_freqs, + H_int, + GAMMA, + GAMMA_PHI, + ) + return pool.submit( + run_simulation_double_pulse_phase, + freq_star, + AMP_DRIVE, + AMP_DRIVE, + PULSE_NS, + FIXED_GAP_NS, + PULSE_NS, + p, + tlist, + init_freqs, + H_int, + GAMMA, + GAMMA_PHI, + ) n_workers = max(1, min(workers, len(param_list))) with ProcessPoolExecutor(n_workers) as pool: @@ -362,8 +421,10 @@ def submit(pool, p): # METADATA # ───────────────────────────────────────────────────────────────────────────── -def build_metadata(freq_star, init_freqs, n_freq, n_gap, n_phase, n_t, - smoke_test=False) -> dict: + +def build_metadata( + freq_star, init_freqs, n_freq, n_gap, n_phase, n_t, smoke_test=False +) -> dict: return { "experiment_id": "experiment_9", "title": "Separation & Relative-Phase Coherent Control (N=2 TLS simulation)", @@ -420,18 +481,33 @@ def build_metadata(freq_star, init_freqs, n_freq, n_gap, n_phase, n_t, }, }, "sweep": { - "panels": {"n": 3, "names": PANEL_NAMES, - "note": "panel_id 0=ring-down map, 1=gap sweep, " - "2=phase sweep."}, - "drive_frequency": {"panel": 0, "n": n_freq, "unit": "GHz", - "range": [F_MIN, F_MAX], - "note": "Swept in panel 0; FREQ_STAR in panels 1,2."}, - "gap": {"panel": 1, "n": n_gap, "unit": "ns", "range": [GAP_LO, GAP_HI], - "fixed_value_panel_2": FIXED_GAP_NS, - "note": "Inter-pulse gap; swept in panel 1, fixed in panel 2."}, - "phase2": {"panel": 2, "n": n_phase, "unit": "rad", - "range": [PHASE_LO, PHASE_HI], - "note": "Relative phase of pulse 2; swept in panel 2, else 0."}, + "panels": { + "n": 3, + "names": PANEL_NAMES, + "note": "panel_id 0=ring-down map, 1=gap sweep, " "2=phase sweep.", + }, + "drive_frequency": { + "panel": 0, + "n": n_freq, + "unit": "GHz", + "range": [F_MIN, F_MAX], + "note": "Swept in panel 0; FREQ_STAR in panels 1,2.", + }, + "gap": { + "panel": 1, + "n": n_gap, + "unit": "ns", + "range": [GAP_LO, GAP_HI], + "fixed_value_panel_2": FIXED_GAP_NS, + "note": "Inter-pulse gap; swept in panel 1, fixed in panel 2.", + }, + "phase2": { + "panel": 2, + "n": n_phase, + "unit": "rad", + "range": [PHASE_LO, PHASE_HI], + "note": "Relative phase of pulse 2; swept in panel 2, else 0.", + }, "amplitude": {"value": AMP_DRIVE, "note": "Both pulses; not swept."}, "time_samples_per_trace": n_t, }, @@ -446,7 +522,7 @@ def build_metadata(freq_star, init_freqs, n_freq, n_gap, n_phase, n_t, "provenance": { "physics_module": "hamiltonian_generator.py (QuTiP-5 API shim only)", "attribution": "Uses code and simulation methods from the " - "Fitzpatrick Lab, Dartmouth College.", + "Fitzpatrick Lab, Dartmouth College.", "no_new_physics": True, "smoke_test": smoke_test, }, @@ -466,58 +542,112 @@ def build_column_info(n_freq, n_gap, n_phase, n_t) -> dict: "run live once in metadata.json, not repeated per row." ), "columns": [ - {"name": "panel_id", "dtype": "float32", "unit": "categorical (0/1/2)", - "role": "coordinate (regime selector)", "measured": False, - "values": [0, 1, 2], - "value_labels": {0: PANEL_NAMES[0], 1: PANEL_NAMES[1], 2: PANEL_NAMES[2]}, - "n_unique": 3, - "description": "0 ring-down map (single pulse, drive freq swept), " - "1 gap sweep (two pulses, phase2=0, gap swept, at FREQ_STAR), " - "2 phase sweep (two pulses, fixed gap, phase2 swept, at FREQ_STAR)."}, - {"name": "drive_frequency_GHz", "dtype": "float32", - "unit": "GHz (== ns^-1, no 2*pi)", - "role": "coordinate (swept input in panel 0)", "measured": False, - "range": [F_MIN, F_MAX], "n_unique_panel_0": n_freq, - "description": "Microwave drive frequency. Swept across 2.0-5.0 GHz in " - "panel 0; held at FREQ_STAR in panels 1 and 2."}, - {"name": "amp1", "dtype": "float32", - "unit": "dimensionless drive amplitude", - "role": "coordinate (fixed)", "measured": False, "value": AMP_DRIVE, - "description": "Amplitude of the first (or only) drive pulse. Constant " - "at 0.12 throughout."}, - {"name": "amp2", "dtype": "float32", - "unit": "dimensionless drive amplitude", - "role": "coordinate", "measured": False, - "description": "Amplitude of the second drive pulse (0.12 in panels 1,2; " - "0 in panel 0 where there is no second pulse)."}, - {"name": "pulse1_ns", "dtype": "float32", "unit": "ns", - "role": "coordinate (pulse timing)", "measured": False, - "description": "Duration of the first pulse. Constant (200 ns)."}, - {"name": "gap_ns", "dtype": "float32", "unit": "ns", - "role": "coordinate (swept input in panel 1)", "measured": False, - "range": [GAP_LO, GAP_HI], - "description": "Idle gap after pulse 1 before pulse 2. Swept 0-800 ns in " - "panel 1; fixed at 150 ns in panel 2; 0 in panel 0 (single pulse)."}, - {"name": "pulse2_ns", "dtype": "float32", "unit": "ns", - "role": "coordinate (pulse timing)", "measured": False, - "description": "Duration of the second pulse (200 ns in panels 1,2; " - "0 in panel 0)."}, - {"name": "phase2_rad", "dtype": "float32", "unit": "rad", - "role": "coordinate (swept input in panel 2)", "measured": False, - "range": [PHASE_LO, PHASE_HI], - "description": "Relative phase of the second pulse. Swept 0-2*pi in " - "panel 2; 0 elsewhere. Decides whether the second pulse's response " - "adds to or cancels the ring-down of the first."}, - {"name": "timestamp_ns", "dtype": "float32", "unit": "ns", - "role": "coordinate (time axis)", "measured": False, - "range": [0.0, float(T_MAX_NS)], "n_unique": n_t, - "description": "Simulation time within the trace. Identical grid for " - "every trace, so traces are directly comparable."}, - {"name": "population", "dtype": "float32", "unit": "dimensionless", - "role": "simulated observable (output)", "measured": False, - "computed_from": "QuTiP mesolve expectation of Sm.dag()*Sm", - "description": "Collective excitation of the N=2 TLS " - "pair, computed by the Lindblad solver. Non-negative. The sole output."}, + { + "name": "panel_id", + "dtype": "float32", + "unit": "categorical (0/1/2)", + "role": "coordinate (regime selector)", + "measured": False, + "values": [0, 1, 2], + "value_labels": { + 0: PANEL_NAMES[0], + 1: PANEL_NAMES[1], + 2: PANEL_NAMES[2], + }, + "n_unique": 3, + "description": "0 ring-down map (single pulse, drive freq swept), " + "1 gap sweep (two pulses, phase2=0, gap swept, at FREQ_STAR), " + "2 phase sweep (two pulses, fixed gap, phase2 swept, at FREQ_STAR).", + }, + { + "name": "drive_frequency_GHz", + "dtype": "float32", + "unit": "GHz (== ns^-1, no 2*pi)", + "role": "coordinate (swept input in panel 0)", + "measured": False, + "range": [F_MIN, F_MAX], + "n_unique_panel_0": n_freq, + "description": "Microwave drive frequency. Swept across 2.0-5.0 GHz in " + "panel 0; held at FREQ_STAR in panels 1 and 2.", + }, + { + "name": "amp1", + "dtype": "float32", + "unit": "dimensionless drive amplitude", + "role": "coordinate (fixed)", + "measured": False, + "value": AMP_DRIVE, + "description": "Amplitude of the first (or only) drive pulse. Constant " + "at 0.12 throughout.", + }, + { + "name": "amp2", + "dtype": "float32", + "unit": "dimensionless drive amplitude", + "role": "coordinate", + "measured": False, + "description": "Amplitude of the second drive pulse (0.12 in panels 1,2; " + "0 in panel 0 where there is no second pulse).", + }, + { + "name": "pulse1_ns", + "dtype": "float32", + "unit": "ns", + "role": "coordinate (pulse timing)", + "measured": False, + "description": "Duration of the first pulse. Constant (200 ns).", + }, + { + "name": "gap_ns", + "dtype": "float32", + "unit": "ns", + "role": "coordinate (swept input in panel 1)", + "measured": False, + "range": [GAP_LO, GAP_HI], + "description": "Idle gap after pulse 1 before pulse 2. Swept 0-800 ns in " + "panel 1; fixed at 150 ns in panel 2; 0 in panel 0 (single pulse).", + }, + { + "name": "pulse2_ns", + "dtype": "float32", + "unit": "ns", + "role": "coordinate (pulse timing)", + "measured": False, + "description": "Duration of the second pulse (200 ns in panels 1,2; " + "0 in panel 0).", + }, + { + "name": "phase2_rad", + "dtype": "float32", + "unit": "rad", + "role": "coordinate (swept input in panel 2)", + "measured": False, + "range": [PHASE_LO, PHASE_HI], + "description": "Relative phase of the second pulse. Swept 0-2*pi in " + "panel 2; 0 elsewhere. Decides whether the second pulse's response " + "adds to or cancels the ring-down of the first.", + }, + { + "name": "timestamp_ns", + "dtype": "float32", + "unit": "ns", + "role": "coordinate (time axis)", + "measured": False, + "range": [0.0, float(T_MAX_NS)], + "n_unique": n_t, + "description": "Simulation time within the trace. Identical grid for " + "every trace, so traces are directly comparable.", + }, + { + "name": "population", + "dtype": "float32", + "unit": "dimensionless", + "role": "simulated observable (output)", + "measured": False, + "computed_from": "QuTiP mesolve expectation of Sm.dag()*Sm", + "description": "Collective excitation of the N=2 TLS " + "pair, computed by the Lindblad solver. Non-negative. The sole output.", + }, ], "computable_not_stored": { "note": "population is the only observable stored here.", @@ -529,6 +659,7 @@ def build_column_info(n_freq, n_gap, n_phase, n_t) -> dict: # I/O # ───────────────────────────────────────────────────────────────────────────── + class _Tee: """Duplicate everything written to stdout into a log file as well, so each run leaves a persistent, human-readable record for reverification (the run @@ -566,7 +697,7 @@ def write_csv(data: np.ndarray, columns: list, path: str, chunk_rows: int = 500_ first = True with open(path, "w", newline="", encoding="utf-8") as fh: for start in tqdm(range(0, n, chunk_rows), desc="CSV", unit="chunk"): - df = pd.DataFrame(data[start: start + chunk_rows], columns=columns) + df = pd.DataFrame(data[start : start + chunk_rows], columns=columns) df.to_csv(fh, index=False, header=first) first = False @@ -575,6 +706,7 @@ def write_csv(data: np.ndarray, columns: list, path: str, chunk_rows: int = 500_ # MAIN # ───────────────────────────────────────────────────────────────────────────── + def build_dataset(output_dir, workers=4, write_csv_file=True, smoke_test=False): """Public entry point. Tees all console output to a run log for reverification, then delegates to _build_dataset_impl.""" @@ -586,12 +718,15 @@ def build_dataset(output_dir, workers=4, write_csv_file=True, smoke_test=False): with open(log_path, "w", encoding="utf-8") as log_fh: sys.stdout = _Tee(original_stdout, log_fh) try: - print(f"# Experiment 9 run log — " - f"{datetime.datetime.now():%Y-%m-%d %H:%M:%S} " - f"({'SMOKE TEST' if smoke_test else 'full run'}, " - f"workers={workers})") - result = _build_dataset_impl(output_dir, workers, write_csv_file, - smoke_test) + print( + f"# Experiment 9 run log — " + f"{datetime.datetime.now():%Y-%m-%d %H:%M:%S} " + f"({'SMOKE TEST' if smoke_test else 'full run'}, " + f"workers={workers})" + ) + result = _build_dataset_impl( + output_dir, workers, write_csv_file, smoke_test + ) finally: sys.stdout = original_stdout @@ -612,17 +747,17 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F else: n_freq, n_gap, n_phase, n_t = N_FREQ, N_GAP, N_PHASE, N_T - freq_axis = np.linspace(F_MIN, F_MAX, n_freq) - gap_sweep = np.linspace(GAP_LO, GAP_HI, n_gap) + freq_axis = np.linspace(F_MIN, F_MAX, n_freq) + gap_sweep = np.linspace(GAP_LO, GAP_HI, n_gap) phase_sweep = np.linspace(PHASE_LO, PHASE_HI, n_phase) - tlist = np.linspace(0.0, T_MAX_NS, n_t) + tlist = np.linspace(0.0, T_MAX_NS, n_t) n_total_rows = (n_freq + n_gap + n_phase) * n_t print(f"\nExperiment 9 dataset builder {'(SMOKE TEST)' if smoke_test else ''}") - print(f"Format : long (one row per time sample)") + print("Format : long (one row per time sample)") print(f"System : N={N_TLS} coupled TLS, Lindblad master equation") - print(f"Panels : 0=ring-down map, 1=gap sweep, 2=phase sweep") + print("Panels : 0=ring-down map, 1=gap sweep, 2=phase sweep") print(f"Drive freqs : {n_freq} ({F_MIN}-{F_MAX} GHz, panel 0)") print(f"Gaps : {n_gap} ({GAP_LO}-{GAP_HI} ns, panel 1)") print(f"Phases : {n_phase} (0-2*pi rad, panel 2)") @@ -635,8 +770,9 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F # ── panel 0: ring-down map + FREQ_STAR ─────────────────────────────────── print("[panel 0] ring-down map ...") - pop_ring = _run_sweep("freq", freq_axis, None, tlist, init_freqs, H_int, - workers, smoke_test) + pop_ring = _run_sweep( + "freq", freq_axis, None, tlist, init_freqs, H_int, workers, smoke_test + ) mask_tail = tlist > PULSE_NS if mask_tail.any(): freq_star = float(freq_axis[np.argmax(pop_ring[:, mask_tail].sum(axis=1))]) @@ -646,13 +782,15 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F # ── panel 1: gap sweep (phase2 = 0) ────────────────────────────────────── print("\n[panel 1] gap sweep ...") - pop_gap = _run_sweep("gap", gap_sweep, freq_star, tlist, init_freqs, H_int, - workers, smoke_test) + pop_gap = _run_sweep( + "gap", gap_sweep, freq_star, tlist, init_freqs, H_int, workers, smoke_test + ) # ── panel 2: phase sweep (fixed gap) ───────────────────────────────────── print("\n[panel 2] phase sweep ...") - pop_phase = _run_sweep("phase", phase_sweep, freq_star, tlist, init_freqs, H_int, - workers, smoke_test) + pop_phase = _run_sweep( + "phase", phase_sweep, freq_star, tlist, init_freqs, H_int, workers, smoke_test + ) # ── assemble the long table ────────────────────────────────────────────── data_all = np.empty((n_total_rows, len(COLUMN_NAMES)), dtype=np.float32) @@ -664,15 +802,16 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, phase2, pop_matrix, param_it end = row_idx + n_t data_all[row_idx:end, COL["panel_id"]] = panel_id data_all[row_idx:end, COL["drive_frequency_GHz"]] = ( - drive_freq[row_i] if np.ndim(drive_freq) else drive_freq) + drive_freq[row_i] if np.ndim(drive_freq) else drive_freq + ) data_all[row_idx:end, COL["amp1"]] = a1 data_all[row_idx:end, COL["amp2"]] = a2 data_all[row_idx:end, COL["pulse1_ns"]] = p1 - data_all[row_idx:end, COL["gap_ns"]] = ( - gap[row_i] if np.ndim(gap) else gap) + data_all[row_idx:end, COL["gap_ns"]] = gap[row_i] if np.ndim(gap) else gap data_all[row_idx:end, COL["pulse2_ns"]] = p2 data_all[row_idx:end, COL["phase2_rad"]] = ( - phase2[row_i] if np.ndim(phase2) else phase2) + phase2[row_i] if np.ndim(phase2) else phase2 + ) data_all[row_idx:end, COL["timestamp_ns"]] = tlist data_all[row_idx:end, COL["population"]] = pop_matrix[row_i] row_idx = end @@ -680,18 +819,39 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, phase2, pop_matrix, param_it # panel 0: single pulse -> amp2=0, gap=0, pulse2=0, phase2=0 emit(0, freq_axis, AMP_DRIVE, 0.0, PULSE_NS, 0.0, 0.0, 0.0, pop_ring, freq_axis) # panel 1: two pulses, phase2=0, gap swept - emit(1, freq_star, AMP_DRIVE, AMP_DRIVE, PULSE_NS, gap_sweep, PULSE_NS, 0.0, - pop_gap, gap_sweep) + emit( + 1, + freq_star, + AMP_DRIVE, + AMP_DRIVE, + PULSE_NS, + gap_sweep, + PULSE_NS, + 0.0, + pop_gap, + gap_sweep, + ) # panel 2: two pulses, gap fixed, phase2 swept - emit(2, freq_star, AMP_DRIVE, AMP_DRIVE, PULSE_NS, FIXED_GAP_NS, PULSE_NS, - phase_sweep, pop_phase, phase_sweep) + emit( + 2, + freq_star, + AMP_DRIVE, + AMP_DRIVE, + PULSE_NS, + FIXED_GAP_NS, + PULSE_NS, + phase_sweep, + pop_phase, + phase_sweep, + ) data_final = data_all[:row_idx] # ── save pickle + JSON sidecars ────────────────────────────────────────── print(f"\nSaving to {pkl_path} ...") - metadata = build_metadata(freq_star, init_freqs, n_freq, n_gap, n_phase, n_t, - smoke_test) + metadata = build_metadata( + freq_star, init_freqs, n_freq, n_gap, n_phase, n_t, smoke_test + ) column_doc = build_column_info(n_freq, n_gap, n_phase, n_t) payload = { @@ -709,14 +869,14 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, phase2, pop_matrix, param_it with open(cols_path, "w", encoding="utf-8") as fh: json.dump(column_doc, fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)})") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") # ── sanity checks ──────────────────────────────────────────────────────── - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] if row_idx != n_total_rows: @@ -729,7 +889,7 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, phase2, pop_matrix, param_it if n_nan or n_inf: errors.append(f" [FAIL] {n_nan} NaN and {n_inf} Inf values present") else: - print(f" [PASS] No NaN or Inf values") + print(" [PASS] No NaN or Inf values") pop = data_final[:, COL["population"]] if pop.min() < -1e-3: @@ -754,16 +914,21 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, phase2, pop_matrix, param_it p2 = data_final[:, COL["panel_id"]] == 2 # Structural invariants that must hold across the three panels. - ok_single = (np.all(data_final[p0, COL["amp2"]] == 0.0) and - np.all(data_final[p0, COL["pulse2_ns"]] == 0.0) and - np.all(data_final[p0, COL["phase2_rad"]] == 0.0)) + ok_single = ( + np.all(data_final[p0, COL["amp2"]] == 0.0) + and np.all(data_final[p0, COL["pulse2_ns"]] == 0.0) + and np.all(data_final[p0, COL["phase2_rad"]] == 0.0) + ) if not ok_single: - errors.append(" [FAIL] panel 0 should be single-pulse (amp2/pulse2/phase2 = 0)") + errors.append( + " [FAIL] panel 0 should be single-pulse (amp2/pulse2/phase2 = 0)" + ) else: print(" [PASS] panel 0 is single-pulse (amp2=pulse2=phase2=0)") - ok_two = (np.all(data_final[p1 | p2, COL["amp2"]] == AMP_DRIVE) and - np.all(data_final[p1 | p2, COL["pulse2_ns"]] == PULSE_NS)) + ok_two = np.all(data_final[p1 | p2, COL["amp2"]] == AMP_DRIVE) and np.all( + data_final[p1 | p2, COL["pulse2_ns"]] == PULSE_NS + ) if not ok_two: errors.append(" [FAIL] panels 1,2 should be two-pulse (amp2, pulse2_ns set)") else: @@ -791,21 +956,23 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, phase2, pop_matrix, param_it if not np.array_equal(reloaded["data"], data_final): errors.append(" [FAIL] Pickle round-trip: reloaded data differs") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly") + print(" [PASS] Pickle round-trip: reloaded array matches exactly") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) return pkl_path @@ -818,14 +985,26 @@ def emit(panel_id, drive_freq, a1, a2, p1, gap, p2, phase2, pop_matrix, param_it description="Build the Experiment 9 simulated dataset (long format) as a " "pickle and a CSV (N=2 TLS separation & relative-phase control)." ) - parser.add_argument("--output_dir", type=str, default=_script_dir, - help="Where to write experiment_9/ (default: next to this script).") - parser.add_argument("--workers", type=int, default=4, - help="Parallel solver processes (real run only).") - parser.add_argument("--no_csv", action="store_true", - help="Skip the CSV, write only the pickle.") - parser.add_argument("--smoke-test", action="store_true", - help="Tiny grids + synthetic population; no qutip required.") + parser.add_argument( + "--output_dir", + type=str, + default=_script_dir, + help="Where to write experiment_9/ (default: next to this script).", + ) + parser.add_argument( + "--workers", + type=int, + default=4, + help="Parallel solver processes (real run only).", + ) + parser.add_argument( + "--no_csv", action="store_true", help="Skip the CSV, write only the pickle." + ) + parser.add_argument( + "--smoke-test", + action="store_true", + help="Tiny grids + synthetic population; no qutip required.", + ) args = parser.parse_args() build_dataset( diff --git a/src/dataset_creation_scripts/09_separation_phase_sweep/hamiltonian_generator.py b/src/dataset_creation_scripts/09_separation_phase_sweep/hamiltonian_generator.py index c3bc75f..715a109 100644 --- a/src/dataset_creation_scripts/09_separation_phase_sweep/hamiltonian_generator.py +++ b/src/dataset_creation_scripts/09_separation_phase_sweep/hamiltonian_generator.py @@ -9,13 +9,15 @@ All earlier APIs continue to work. """ + import numpy as np from qutip import qeye, destroy, tensor, mesolve # QuTiP-5 API compat: the QuTiP-4 `Options(progress_bar=None)` object was removed; # solver options are now a plain dict. progress_bar is a display flag only, so # this changes nothing about the physics. -_opts_no_bar = {"progress_bar": None} # silence progress-bar output +_opts_no_bar = {"progress_bar": None} # silence progress-bar output + # ───────── TLS operator factory ───────── def _tls_ops(N): @@ -26,7 +28,7 @@ def _tls_ops(N): sz1 = 2 * sm1.dag() * sm1 - I2 def T(op, k): - return tensor([I2]*k + [op] + [I2]*(N-k-1)) + return tensor([I2] * k + [op] + [I2] * (N - k - 1)) sm = [T(sm1, k) for k in range(N)] sx = [T(sx1, k) for k in range(N)] @@ -36,35 +38,49 @@ def T(op, k): # ───────── single-pulse (pop, Sp, Sm) ───────── -def run_simulation_single_pulse_full(freq, amp, tlist, init_freqs, interactions, - gamma, gamma_phi, pulse_ns, T_ramp_ns=1.0): +def run_simulation_single_pulse_full( + freq, + amp, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + pulse_ns, + T_ramp_ns=1.0, +): N = len(init_freqs) sm, sx, _, sz = _tls_ops(N) Sp = sum(s.dag() for s in sm) Sm = sum(sm) Sx = sum(sx) - H_tls = sum(init_freqs[j]*sz[j]/2 for j in range(N)) + H_tls = sum(init_freqs[j] * sz[j] / 2 for j in range(N)) def d_coeff(t, _=None): if t < T_ramp_ns: - env = 0.5*(1-np.cos(np.pi*t/T_ramp_ns)) - return amp*env*np.cos(freq*t) + env = 0.5 * (1 - np.cos(np.pi * t / T_ramp_ns)) + return amp * env * np.cos(freq * t) elif t < pulse_ns: - return amp*np.cos(freq*t) + return amp * np.cos(freq * t) return 0.0 - H = [H_tls + interactions, [Sx, d_coeff]] - ψ0 = (H_tls + interactions).groundstate()[1] + H = [H_tls + interactions, [Sx, d_coeff]] + ψ0 = (H_tls + interactions).groundstate()[1] - c_ops = [np.sqrt(gamma)*Sm] + c_ops = [np.sqrt(gamma) * Sm] if gamma_phi: _, _, _, szs = _tls_ops(N) - c_ops += [np.sqrt(gamma_phi)*z for z in szs] + c_ops += [np.sqrt(gamma_phi) * z for z in szs] pop, sp, sm_exp = mesolve( # QuTiP-5 API compat: e_ops passed by keyword (now keyword-only). - H, ψ0, tlist, c_ops, e_ops=[Sm.dag()*Sm, Sp, Sm], options=_opts_no_bar + H, + ψ0, + tlist, + c_ops, + e_ops=[Sm.dag() * Sm, Sp, Sm], + options=_opts_no_bar, ).expect return np.real_if_close(pop), sp, sm_exp @@ -75,25 +91,35 @@ def run_simulation_single_pulse(*args, **kw): # ───────── two-pulse with phase shift φ₂ ───────── -def run_simulation_double_pulse_phase(freq, amp1, amp2, - pulse1_ns, gap_ns, pulse2_ns, - phase2_rad, - tlist, init_freqs, interactions, - gamma, gamma_phi, T_ramp_ns=1.0): +def run_simulation_double_pulse_phase( + freq, + amp1, + amp2, + pulse1_ns, + gap_ns, + pulse2_ns, + phase2_rad, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + T_ramp_ns=1.0, +): N = len(init_freqs) sm, sx, _, sz = _tls_ops(N) Sx, Sm = sum(sx), sum(sm) - H_tls = sum(init_freqs[j]*sz[j]/2 for j in range(N)) + H_tls = sum(init_freqs[j] * sz[j] / 2 for j in range(N)) t_start_p2 = pulse1_ns + gap_ns def d_coeff(t, _=None): # pulse-1 if t < T_ramp_ns: - env = 0.5*(1-np.cos(np.pi*t/T_ramp_ns)) - return amp1*env*np.cos(freq*t) + env = 0.5 * (1 - np.cos(np.pi * t / T_ramp_ns)) + return amp1 * env * np.cos(freq * t) elif t < pulse1_ns: - return amp1*np.cos(freq*t) + return amp1 * np.cos(freq * t) # gap if t < t_start_p2: @@ -102,48 +128,77 @@ def d_coeff(t, _=None): # pulse-2 (phase-shifted) t2 = t - t_start_p2 if t2 < T_ramp_ns: - env = 0.5*(1-np.cos(np.pi*t2/T_ramp_ns)) - return amp2*env*np.cos(freq*t + phase2_rad) + env = 0.5 * (1 - np.cos(np.pi * t2 / T_ramp_ns)) + return amp2 * env * np.cos(freq * t + phase2_rad) elif t2 < pulse2_ns: - return amp2*np.cos(freq*t + phase2_rad) + return amp2 * np.cos(freq * t + phase2_rad) return 0.0 - H = [H_tls + interactions, [Sx, d_coeff]] - ψ0 = (H_tls + interactions).groundstate()[1] - c_ops = [np.sqrt(gamma)*Sm] + H = [H_tls + interactions, [Sx, d_coeff]] + ψ0 = (H_tls + interactions).groundstate()[1] + c_ops = [np.sqrt(gamma) * Sm] if gamma_phi: _, _, _, szs = _tls_ops(N) - c_ops += [np.sqrt(gamma_phi)*z for z in szs] + c_ops += [np.sqrt(gamma_phi) * z for z in szs] pop = mesolve( # QuTiP-5 API compat: e_ops passed by keyword (now keyword-only). - H, ψ0, tlist, c_ops, e_ops=[Sm.dag()*Sm], options=_opts_no_bar + H, + ψ0, + tlist, + c_ops, + e_ops=[Sm.dag() * Sm], + options=_opts_no_bar, ).expect[0] return np.real_if_close(pop) # backwards-compat API (phase=0) -def run_simulation_double_pulse(freq, amp1, amp2, - pulse1_ns, gap_ns, pulse2_ns, - tlist, init_freqs, interactions, - gamma, gamma_phi, T_ramp_ns=1.0): +def run_simulation_double_pulse( + freq, + amp1, + amp2, + pulse1_ns, + gap_ns, + pulse2_ns, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + T_ramp_ns=1.0, +): return run_simulation_double_pulse_phase( - freq, amp1, amp2, pulse1_ns, gap_ns, pulse2_ns, 0.0, - tlist, init_freqs, interactions, gamma, gamma_phi, T_ramp_ns + freq, + amp1, + amp2, + pulse1_ns, + gap_ns, + pulse2_ns, + 0.0, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + T_ramp_ns, ) # ───────── random couplings helper ───────── -def build_spin_spin_interactions_random_distribution(N, J_min, J_max, - alpha_x=1.0, alpha_y=0.0, alpha_z=0.0): +def build_spin_spin_interactions_random_distribution( + N, J_min, J_max, alpha_x=1.0, alpha_y=0.0, alpha_z=0.0 +): _, sx, sy, sz = _tls_ops(N) J = np.random.uniform(J_min, J_max, size=(N, N)) J = np.tril(J, -1) + np.tril(J, -1).T H_int = 0 for i in range(N): - for j in range(i+1, N): + for j in range(i + 1, N): Jij = J[i, j] - H_int += Jij*(alpha_x*sx[i]*sx[j] + - alpha_y*sy[i]*sy[j] + - alpha_z*sz[i]*sz[j]) + H_int += Jij * ( + alpha_x * sx[i] * sx[j] + + alpha_y * sy[i] * sy[j] + + alpha_z * sz[i] * sz[j] + ) return H_int diff --git a/src/dataset_creation_scripts/10_sharpening_physics/experiment_10_dataset_creation.py b/src/dataset_creation_scripts/10_sharpening_physics/experiment_10_dataset_creation.py index 740393a..f2b285f 100644 --- a/src/dataset_creation_scripts/10_sharpening_physics/experiment_10_dataset_creation.py +++ b/src/dataset_creation_scripts/10_sharpening_physics/experiment_10_dataset_creation.py @@ -135,33 +135,40 @@ # Drive-frequency sweep. Treated numerically as GHz == ns^-1. F_MIN, F_MAX = 3.0, 5.0 -N_FREQ = 400 +N_FREQ = 400 # Time grid: 0 to 1000 ns, 1000 points. T_MAX_NS, N_T = 1000, 1000 # The N = 2 TLS. Frequencies hard-coded; coupling is built (J=0.05 XX) but the # source solver DROPS it -> the TLS are effectively uncoupled (see header). -N_TLS = 2 +N_TLS = 2 INIT_FREQS = [4.0, 4.1] -J_COUP = 0.05 +J_COUP = 0.05 # Dissipation and drive amplitude. GAMMA, GAMMA_PHI = 0.002, 0.0 -DRIVE_AMPL = 0.1 +DRIVE_AMPL = 0.1 # Phase-V FFT settings. -FFT_VIEW_MHZ = 150.0 +FFT_VIEW_MHZ = 150.0 PHASEV_WIN_START_NS = 0.0 -PHASEV_WIN_STOP_NS = 400.0 +PHASEV_WIN_STOP_NS = 400.0 # ── table schemas ──────────────────────────────────────────────────────────── -POSTPULSE_COLUMNS = ["pulse_duration_ns", "drive_frequency_GHz", - "postpulse_population"] -POPTRACE_COLUMNS = ["pulse_duration_ns", "drive_frequency_GHz", - "timestamp_ns", "population"] -PHASEV_COLUMNS = ["pulse_duration_ns", "drive_frequency_GHz", - "fft_frequency_MHz", "normalized_phase_fft"] +POSTPULSE_COLUMNS = ["pulse_duration_ns", "drive_frequency_GHz", "postpulse_population"] +POPTRACE_COLUMNS = [ + "pulse_duration_ns", + "drive_frequency_GHz", + "timestamp_ns", + "population", +] +PHASEV_COLUMNS = [ + "pulse_duration_ns", + "drive_frequency_GHz", + "fft_frequency_MHz", + "normalized_phase_fft", +] PICKLE_PROTOCOL = 4 @@ -175,12 +182,14 @@ def top_pulse_indices(pulses): # ENSEMBLE # ───────────────────────────────────────────────────────────────────────────── + def build_H_int(): """Build the (nominal) XX coupling operator. NOTE: the solver discards this operator (see the UNCOUPLED note in the module docstring); it is built here only so the real solver path receives the same arguments as the original implementation.""" from hamiltonian_generator import build_spin_spin_interactions_random_distribution + return build_spin_spin_interactions_random_distribution( N_TLS, J_COUP, J_COUP, alpha_x=1.0, alpha_y=0.0, alpha_z=0.0 ) @@ -190,14 +199,17 @@ def build_H_int(): # SYNTHETIC ARRAYS — smoke-test stand-in for the QuTiP solver (no physics). # ───────────────────────────────────────────────────────────────────────────── + def _synthetic_full(freq, amp, tlist, init_freqs, pulse_ns): """Return (sp1, sp2, pop) surrogates. NOT physical; plumbing only.""" detune = np.min(np.abs(np.asarray(init_freqs) - freq)) lorentz = 1.0 / (1.0 + (detune / 0.15) ** 2) - drive = amp * lorentz * np.sin( - np.pi * np.clip(tlist / max(pulse_ns, 1e-9), 0, 1)) ** 2 - tail = np.where(tlist > pulse_ns, - amp * lorentz * np.exp(-GAMMA * (tlist - pulse_ns)), 0.0) + drive = ( + amp * lorentz * np.sin(np.pi * np.clip(tlist / max(pulse_ns, 1e-9), 0, 1)) ** 2 + ) + tail = np.where( + tlist > pulse_ns, amp * lorentz * np.exp(-GAMMA * (tlist - pulse_ns)), 0.0 + ) pop = np.clip(np.where(tlist <= pulse_ns, drive, tail), 0.0, None) env = np.exp(-GAMMA * tlist) sp1 = amp * env * np.exp(1j * init_freqs[0] * tlist * 0.05) @@ -217,16 +229,32 @@ def _init_worker(tlist, init_freqs, H_int, gamma, gamma_phi, amp): so each task only ships (freq, pulse_ns, j_post, want_full).""" global _G from hamiltonian_generator import run_simulation_for_frequency - _G = dict(tlist=tlist, init_freqs=init_freqs, H_int=H_int, gamma=gamma, - gamma_phi=gamma_phi, amp=amp, fn=run_simulation_for_frequency) + + _G = dict( + tlist=tlist, + init_freqs=init_freqs, + H_int=H_int, + gamma=gamma, + gamma_phi=gamma_phi, + amp=amp, + fn=run_simulation_for_frequency, + ) def _worker(task): """Return (postpulse_scalar, sp1, sp2, pop) — sp1/sp2/pop only when want_full, to keep inter-process traffic small for the non-representative pulses.""" freq, pulse_ns, j_post, want_full = task - sp1, sp2, pop = _G["fn"](freq, _G["tlist"], _G["init_freqs"], _G["H_int"], - _G["gamma"], _G["gamma_phi"], _G["amp"], pulse_ns) + sp1, sp2, pop = _G["fn"]( + freq, + _G["tlist"], + _G["init_freqs"], + _G["H_int"], + _G["gamma"], + _G["gamma_phi"], + _G["amp"], + pulse_ns, + ) pop = np.asarray(pop, dtype=np.float64).real scalar = float(pop[j_post]) if want_full: @@ -265,8 +293,10 @@ def compute_phaseV(phi_buffer, tlist, dt): # METADATA # ───────────────────────────────────────────────────────────────────────────── -def build_metadata(pulses, top_idx, freq_axis, tlist, fft_freq_MHz, - smoke_test=False) -> dict: + +def build_metadata( + pulses, top_idx, freq_axis, tlist, fft_freq_MHz, smoke_test=False +) -> dict: return { "experiment_id": "experiment_10", "title": "Ring-down & Phase Sharpening vs Pulse Duration (N=2 TLS simulation)", @@ -297,10 +327,10 @@ def build_metadata(pulses, top_idx, freq_axis, tlist, fft_freq_MHz, "observables": { "population": " collective excitation (dimensionless)", "phase_difference": "phi(t) = arg - arg, " - "wrapped to (-pi, pi]", + "wrapped to (-pi, pi]", }, "drive": "square pulse, amp*cos(freq*t) for t <= pulse_ns (no ramp), " - "per the active drive definition", + "per the active drive definition", "unit_convention": ( "Frequencies used numerically as inverse nanoseconds (GHz == ns^-1) " "with NO factor of 2*pi. Amplitude dimensionless; phase in radians; " @@ -317,16 +347,21 @@ def build_metadata(pulses, top_idx, freq_axis, tlist, fft_freq_MHz, }, "grids": { "pulse_durations_ns": { - "n": len(pulses), "min": float(pulses[0]), "max": float(pulses[-1]), + "n": len(pulses), + "min": float(pulses[0]), + "max": float(pulses[-1]), "formula": "linspace(20, 200, 50)", }, "representative_pulses_ns": [round(float(pulses[i]), 4) for i in top_idx], "drive_frequency_GHz": { - "n": len(freq_axis), "range": [F_MIN, F_MAX], + "n": len(freq_axis), + "range": [F_MIN, F_MAX], "formula": "linspace(3.0, 5.0, 400)", }, "time_grid": { - "t_min_ns": 0.0, "t_max_ns": T_MAX_NS, "n_samples": len(tlist), + "t_min_ns": 0.0, + "t_max_ns": T_MAX_NS, + "n_samples": len(tlist), "one_sample_ns": float(tlist[1] - tlist[0]) if len(tlist) > 1 else 0.0, }, "phaseV": { @@ -337,17 +372,23 @@ def build_metadata(pulses, top_idx, freq_axis, tlist, fft_freq_MHz, }, }, "tables": { - "postpulse_map": {"columns": POSTPULSE_COLUMNS, - "n_rows": len(pulses) * len(freq_axis)}, - "population_traces": {"columns": POPTRACE_COLUMNS, - "n_rows": len(top_idx) * len(freq_axis) * len(tlist)}, - "phaseV_fft": {"columns": PHASEV_COLUMNS, - "n_rows": len(top_idx) * len(freq_axis) * len(fft_freq_MHz)}, + "postpulse_map": { + "columns": POSTPULSE_COLUMNS, + "n_rows": len(pulses) * len(freq_axis), + }, + "population_traces": { + "columns": POPTRACE_COLUMNS, + "n_rows": len(top_idx) * len(freq_axis) * len(tlist), + }, + "phaseV_fft": { + "columns": PHASEV_COLUMNS, + "n_rows": len(top_idx) * len(freq_axis) * len(fft_freq_MHz), + }, }, "provenance": { "physics_module": "hamiltonian_generator.py (QuTiP-5 API shim only)", "attribution": "Uses code and simulation methods from the " - "Fitzpatrick Lab, Dartmouth College.", + "Fitzpatrick Lab, Dartmouth College.", "no_new_physics": True, "reproduces_uncoupled_quirk": True, "smoke_test": smoke_test, @@ -365,53 +406,100 @@ def build_column_info(pulses, top_idx, freq_axis, tlist, fft_freq_MHz) -> dict: "format": "long (one row per pulse x drive frequency)", "n_rows": len(pulses) * len(freq_axis), "columns": [ - {"name": "pulse_duration_ns", "dtype": "float32", "unit": "ns", - "role": "coordinate", "description": "Drive pulse duration " - "(20-200 ns, 50 values)."}, - {"name": "drive_frequency_GHz", "dtype": "float32", - "unit": "GHz (== ns^-1)", "role": "coordinate", - "description": "Drive frequency (3.0-5.0 GHz, 400 values)."}, - {"name": "postpulse_population", "dtype": "float32", - "unit": "dimensionless", "role": "observable", - "description": " sampled at the first time " - "sample after the pulse ends (first t > pulse_duration_ns)."}, + { + "name": "pulse_duration_ns", + "dtype": "float32", + "unit": "ns", + "role": "coordinate", + "description": "Drive pulse duration " + "(20-200 ns, 50 values).", + }, + { + "name": "drive_frequency_GHz", + "dtype": "float32", + "unit": "GHz (== ns^-1)", + "role": "coordinate", + "description": "Drive frequency (3.0-5.0 GHz, 400 values).", + }, + { + "name": "postpulse_population", + "dtype": "float32", + "unit": "dimensionless", + "role": "observable", + "description": " sampled at the first time " + "sample after the pulse ends (first t > pulse_duration_ns).", + }, ], }, "population_traces": { "format": "long (one row per representative pulse x freq x time)", "n_rows": len(top_idx) * len(freq_axis) * len(tlist), "columns": [ - {"name": "pulse_duration_ns", "dtype": "float32", "unit": "ns", - "role": "coordinate", "description": "One of the three " - "representative pulses (shortest, middle, longest)."}, - {"name": "drive_frequency_GHz", "dtype": "float32", - "unit": "GHz", "role": "coordinate", - "description": "Drive frequency (3.0-5.0 GHz, 400 values)."}, - {"name": "timestamp_ns", "dtype": "float32", "unit": "ns", - "role": "coordinate (time axis)", - "description": "Simulation time (0-1000 ns, 1000 samples)."}, - {"name": "population", "dtype": "float32", "unit": "dimensionless", - "role": "observable", - "description": "Full time trace."}, + { + "name": "pulse_duration_ns", + "dtype": "float32", + "unit": "ns", + "role": "coordinate", + "description": "One of the three " + "representative pulses (shortest, middle, longest).", + }, + { + "name": "drive_frequency_GHz", + "dtype": "float32", + "unit": "GHz", + "role": "coordinate", + "description": "Drive frequency (3.0-5.0 GHz, 400 values).", + }, + { + "name": "timestamp_ns", + "dtype": "float32", + "unit": "ns", + "role": "coordinate (time axis)", + "description": "Simulation time (0-1000 ns, 1000 samples).", + }, + { + "name": "population", + "dtype": "float32", + "unit": "dimensionless", + "role": "observable", + "description": "Full time trace.", + }, ], }, "phaseV_fft": { "format": "long (one row per representative pulse x freq x fft bin)", "n_rows": len(top_idx) * len(freq_axis) * len(fft_freq_MHz), "columns": [ - {"name": "pulse_duration_ns", "dtype": "float32", "unit": "ns", - "role": "coordinate", "description": "Representative pulse."}, - {"name": "drive_frequency_GHz", "dtype": "float32", - "unit": "GHz", "role": "coordinate", - "description": "Drive frequency (3.0-5.0 GHz, 400 values)."}, - {"name": "fft_frequency_MHz", "dtype": "float32", "unit": "MHz", - "role": "coordinate (FFT axis)", - "description": "Positive FFT frequency of phi(t), 0-150 MHz."}, - {"name": "normalized_phase_fft", "dtype": "float32", - "unit": "dimensionless (0-1)", "role": "observable", - "description": "Row-wise normalized |FFT[phi(t)]| over the " - "0-400 ns window, where phi = arg - arg. " - "Normalized to peak 1 per drive frequency."}, + { + "name": "pulse_duration_ns", + "dtype": "float32", + "unit": "ns", + "role": "coordinate", + "description": "Representative pulse.", + }, + { + "name": "drive_frequency_GHz", + "dtype": "float32", + "unit": "GHz", + "role": "coordinate", + "description": "Drive frequency (3.0-5.0 GHz, 400 values).", + }, + { + "name": "fft_frequency_MHz", + "dtype": "float32", + "unit": "MHz", + "role": "coordinate (FFT axis)", + "description": "Positive FFT frequency of phi(t), 0-150 MHz.", + }, + { + "name": "normalized_phase_fft", + "dtype": "float32", + "unit": "dimensionless (0-1)", + "role": "observable", + "description": "Row-wise normalized |FFT[phi(t)]| over the " + "0-400 ns window, where phi = arg - arg. " + "Normalized to peak 1 per drive frequency.", + }, ], }, }, @@ -422,6 +510,7 @@ def build_column_info(pulses, top_idx, freq_axis, tlist, fft_freq_MHz) -> dict: # I/O # ───────────────────────────────────────────────────────────────────────────── + class _Tee: """Duplicate stdout into a log file for reverification. tqdm uses stderr, so the log stays clean (summary + sanity-check PASS/FAIL lines).""" @@ -451,12 +540,14 @@ def load_pickle(path): def write_csv(data, columns, path, chunk_rows=500_000): import pandas as pd + n = len(data) first = True with open(path, "w", newline="", encoding="utf-8") as fh: - for start in tqdm(range(0, n, chunk_rows), desc=os.path.basename(path), - unit="chunk"): - df = pd.DataFrame(data[start: start + chunk_rows], columns=columns) + for start in tqdm( + range(0, n, chunk_rows), desc=os.path.basename(path), unit="chunk" + ): + df = pd.DataFrame(data[start : start + chunk_rows], columns=columns) df.to_csv(fh, index=False, header=first) first = False @@ -465,11 +556,12 @@ def write_csv(data, columns, path, chunk_rows=500_000): # SOLVE # ───────────────────────────────────────────────────────────────────────────── + def run_all(pulses, top_idx, freq_axis, tlist, workers, smoke_test): """Run the full pulse x frequency sweep, returning: - postpulse_map (n_pulses, n_freq) - top_pop_maps dict pulse_idx -> (n_freq, n_t) - top_phi_maps dict pulse_idx -> (n_freq, n_t) wrapped phase + postpulse_map (n_pulses, n_freq) + top_pop_maps dict pulse_idx -> (n_freq, n_t) + top_phi_maps dict pulse_idx -> (n_freq, n_t) wrapped phase """ n_pulses, n_freq, n_t = len(pulses), len(freq_axis), len(tlist) postpulse_map = np.zeros((n_pulses, n_freq), dtype=np.float32) @@ -489,8 +581,9 @@ def wrap(sp1, sp2): j_post = first_index_after(tlist, float(pulse_ns)) want_full = pi in top_set for fi, f in enumerate(freq_axis): - sp1, sp2, pop = _synthetic_full(f, DRIVE_AMPL, tlist, INIT_FREQS, - float(pulse_ns)) + sp1, sp2, pop = _synthetic_full( + f, DRIVE_AMPL, tlist, INIT_FREQS, float(pulse_ns) + ) postpulse_map[pi, fi] = pop[j_post] if want_full: top_pop_maps[pi][fi] = pop.astype(np.float32) @@ -499,6 +592,7 @@ def wrap(sp1, sp2): # ── real physics: one shared pool, all (pulse, freq) tasks ──────────────── from concurrent.futures import ProcessPoolExecutor, as_completed + H_int = build_H_int() tasks = [] @@ -510,7 +604,8 @@ def wrap(sp1, sp2): n_workers = max(1, min(workers, len(tasks))) with ProcessPoolExecutor( - n_workers, initializer=_init_worker, + n_workers, + initializer=_init_worker, initargs=(tlist, INIT_FREQS, H_int, GAMMA, GAMMA_PHI, DRIVE_AMPL), ) as pool: futures = {pool.submit(_worker, t[2]): (t[0], t[1]) for t in tasks} @@ -528,6 +623,7 @@ def wrap(sp1, sp2): # TABLE ASSEMBLY # ───────────────────────────────────────────────────────────────────────────── + def assemble_postpulse(pulses, freq_axis, postpulse_map): n = len(pulses) * len(freq_axis) data = np.empty((n, 3), dtype=np.float32) @@ -579,6 +675,7 @@ def assemble_phaseV(pulses, top_idx, freq_axis, top_phi_maps, tlist, dt): # MAIN # ───────────────────────────────────────────────────────────────────────────── + def build_dataset(output_dir, workers=4, write_csv_file=True, smoke_test=False): """Public entry point. Tees stdout to a reverification log, then builds.""" exp_dir = os.path.join(output_dir, "experiment_10") @@ -589,12 +686,15 @@ def build_dataset(output_dir, workers=4, write_csv_file=True, smoke_test=False): with open(log_path, "w", encoding="utf-8") as log_fh: sys.stdout = _Tee(original_stdout, log_fh) try: - print(f"# Experiment 10 run log — " - f"{datetime.datetime.now():%Y-%m-%d %H:%M:%S} " - f"({'SMOKE TEST' if smoke_test else 'full run'}, " - f"workers={workers})") - result = _build_dataset_impl(output_dir, workers, write_csv_file, - smoke_test) + print( + f"# Experiment 10 run log — " + f"{datetime.datetime.now():%Y-%m-%d %H:%M:%S} " + f"({'SMOKE TEST' if smoke_test else 'full run'}, " + f"workers={workers})" + ) + result = _build_dataset_impl( + output_dir, workers, write_csv_file, smoke_test + ) finally: sys.stdout = original_stdout @@ -612,39 +712,42 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F else: n_pulses, n_freq, n_t = N_PULSES, N_FREQ, N_T - pulses = np.linspace(PULSE_MIN_NS, PULSE_MAX_NS, n_pulses) + pulses = np.linspace(PULSE_MIN_NS, PULSE_MAX_NS, n_pulses) freq_axis = np.linspace(F_MIN, F_MAX, n_freq) - tlist = np.linspace(0.0, T_MAX_NS, n_t) - dt = float(tlist[1] - tlist[0]) if n_t > 1 else 1.0 - top_idx = top_pulse_indices(pulses) + tlist = np.linspace(0.0, T_MAX_NS, n_t) + dt = float(tlist[1] - tlist[0]) if n_t > 1 else 1.0 + top_idx = top_pulse_indices(pulses) print(f"\nExperiment 10 dataset builder {'(SMOKE TEST)' if smoke_test else ''}") print(f"UNCOUPLED N=2 TLS (source quirk reproduced) at {INIT_FREQS} GHz") print(f"Pulses : {n_pulses} ({PULSE_MIN_NS}-{PULSE_MAX_NS} ns)") - print(f" representative: {[round(float(pulses[i]),2) for i in top_idx]} ns") + print(f" representative: {[round(float(pulses[i]), 2) for i in top_idx]} ns") print(f"Drive freqs : {n_freq} ({F_MIN}-{F_MAX} GHz)") print(f"Time samples : {n_t} (0-{T_MAX_NS} ns)") - print(f"Total solves : {n_pulses} x {n_freq} = {n_pulses*n_freq:,}\n") + print(f"Total solves : {n_pulses} x {n_freq} = {n_pulses * n_freq:,}\n") # ── run the sweep ──────────────────────────────────────────────────────── postpulse_map, top_pop_maps, top_phi_maps = run_all( - pulses, top_idx, freq_axis, tlist, workers, smoke_test) + pulses, top_idx, freq_axis, tlist, workers, smoke_test + ) # ── assemble the three tables ──────────────────────────────────────────── print("\nAssembling tables ...") pp_data = assemble_postpulse(pulses, freq_axis, postpulse_map) pt_data = assemble_poptraces(pulses, top_idx, freq_axis, tlist, top_pop_maps) - pv_data, fft_freq_MHz = assemble_phaseV(pulses, top_idx, freq_axis, - top_phi_maps, tlist, dt) + pv_data, fft_freq_MHz = assemble_phaseV( + pulses, top_idx, freq_axis, top_phi_maps, tlist, dt + ) tables = { - "postpulse_map": (pp_data, POSTPULSE_COLUMNS), + "postpulse_map": (pp_data, POSTPULSE_COLUMNS), "population_traces": (pt_data, POPTRACE_COLUMNS), - "phaseV_fft": (pv_data, PHASEV_COLUMNS), + "phaseV_fft": (pv_data, PHASEV_COLUMNS), } - metadata = build_metadata(pulses, top_idx, freq_axis, tlist, fft_freq_MHz, - smoke_test) + metadata = build_metadata( + pulses, top_idx, freq_axis, tlist, fft_freq_MHz, smoke_test + ) column_doc = build_column_info(pulses, top_idx, freq_axis, tlist, fft_freq_MHz) # ── save pickles + JSON sidecars ───────────────────────────────────────── @@ -652,12 +755,19 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F pkl_paths = {} for name, (data, cols) in tables.items(): pkl_path = os.path.join(output_dir, f"experiment_10_{name}.pkl") - payload = {"data": data, "columns": cols, "table": name, - "attrs": metadata, "column_doc": column_doc["tables"][name]} + payload = { + "data": data, + "columns": cols, + "table": name, + "attrs": metadata, + "column_doc": column_doc["tables"][name], + } save_pickle(payload, pkl_path) pkl_paths[name] = pkl_path - print(f" {name:18s}: {data.shape} -> {pkl_path} " - f"({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print( + f" {name:18s}: {data.shape} -> {pkl_path} " + f"({os.path.getsize(pkl_path) / 1024**2:.1f} MB)" + ) meta_path = os.path.join(script_dir, "metadata.json") cols_path = os.path.join(script_dir, "column_info.json") @@ -680,8 +790,10 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F } for name, (data, cols) in tables.items(): if data.shape[0] != expected[name]: - errors.append(f" [FAIL] {name} rows: got {data.shape[0]:,}, " - f"expected {expected[name]:,}") + errors.append( + f" [FAIL] {name} rows: got {data.shape[0]:,}, " + f"expected {expected[name]:,}" + ) elif int(np.isnan(data).sum()) or int(np.isinf(data).sum()): errors.append(f" [FAIL] {name} contains NaN/Inf") else: @@ -691,23 +803,31 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F if pop_col.min() < -1e-3: errors.append(f" [FAIL] population negative: {pop_col.min():.4f}") else: - print(f" [PASS] population non-negative: " - f"[{pop_col.min():.4f}, {pop_col.max():.4f}]") + print( + f" [PASS] population non-negative: " + f"[{pop_col.min():.4f}, {pop_col.max():.4f}]" + ) pp_col = pp_data[:, POSTPULSE_COLUMNS.index("postpulse_population")] if pp_col.min() < -1e-3: errors.append(f" [FAIL] postpulse_population negative: {pp_col.min():.4f}") else: - print(f" [PASS] postpulse_population non-negative: " - f"[{pp_col.min():.4f}, {pp_col.max():.4f}]") + print( + f" [PASS] postpulse_population non-negative: " + f"[{pp_col.min():.4f}, {pp_col.max():.4f}]" + ) pv_col = pv_data[:, PHASEV_COLUMNS.index("normalized_phase_fft")] if pv_col.min() < -1e-6 or pv_col.max() > 1.0 + 1e-4: - errors.append(f" [FAIL] normalized_phase_fft outside [0,1]: " - f"[{pv_col.min():.4f}, {pv_col.max():.4f}]") + errors.append( + f" [FAIL] normalized_phase_fft outside [0,1]: " + f"[{pv_col.min():.4f}, {pv_col.max():.4f}]" + ) else: - print(f" [PASS] normalized_phase_fft in [0,1]: " - f"[{pv_col.min():.4f}, {pv_col.max():.4f}]") + print( + f" [PASS] normalized_phase_fft in [0,1]: " + f"[{pv_col.min():.4f}, {pv_col.max():.4f}]" + ) if len(np.unique(pp_data[:, 0])) != n_pulses: errors.append(" [FAIL] postpulse_map pulse count mismatch") @@ -734,8 +854,10 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F for name, (data, cols) in tables.items(): csv_path = os.path.join(output_dir, f"experiment_10_{name}.csv") write_csv(data, cols, csv_path) - print(f" {name:18s}: {os.path.getsize(csv_path)/1024**2:.1f} MB, " - f"{data.shape[0]:,} rows") + print( + f" {name:18s}: {os.path.getsize(csv_path) / 1024**2:.1f} MB, " + f"{data.shape[0]:,} rows" + ) return pkl_paths @@ -747,14 +869,26 @@ def _build_dataset_impl(output_dir, workers=4, write_csv_file=True, smoke_test=F "tables; N=2 TLS ring-down & phase sharpening, including the uncoupled " "quirk)." ) - parser.add_argument("--output_dir", type=str, default=_script_dir, - help="Where to write experiment_10/ (default: next to this script).") - parser.add_argument("--workers", type=int, default=4, - help="Parallel solver processes (real run only).") - parser.add_argument("--no_csv", action="store_true", - help="Skip the CSVs, write only the pickles.") - parser.add_argument("--smoke-test", action="store_true", - help="Tiny grids + synthetic arrays; no qutip required.") + parser.add_argument( + "--output_dir", + type=str, + default=_script_dir, + help="Where to write experiment_10/ (default: next to this script).", + ) + parser.add_argument( + "--workers", + type=int, + default=4, + help="Parallel solver processes (real run only).", + ) + parser.add_argument( + "--no_csv", action="store_true", help="Skip the CSVs, write only the pickles." + ) + parser.add_argument( + "--smoke-test", + action="store_true", + help="Tiny grids + synthetic arrays; no qutip required.", + ) args = parser.parse_args() build_dataset( diff --git a/src/dataset_creation_scripts/10_sharpening_physics/hamiltonian_generator.py b/src/dataset_creation_scripts/10_sharpening_physics/hamiltonian_generator.py index 4f3c92c..7528a4d 100644 --- a/src/dataset_creation_scripts/10_sharpening_physics/hamiltonian_generator.py +++ b/src/dataset_creation_scripts/10_sharpening_physics/hamiltonian_generator.py @@ -13,8 +13,14 @@ """ import numpy as np -from qutip import qeye, destroy, tensor, mesolve -from qutip import * +from qutip import ( + qeye, + destroy, + tensor, + mesolve, + sigmap, +) + # ──────────────────── TLS operators ──────────────────── def _tls_ops(N): @@ -24,8 +30,8 @@ def _tls_ops(N): sy1 = -1j * (sm1 - sm1.dag()) sz1 = 2 * sm1.dag() * sm1 - I2 - def T(op, k): # tensor helper - return tensor([I2]*k + [op] + [I2]*(N-k-1)) + def T(op, k): # tensor helper + return tensor([I2] * k + [op] + [I2] * (N - k - 1)) sm = [T(sm1, k) for k in range(N)] sx = [T(sx1, k) for k in range(N)] @@ -35,26 +41,38 @@ def T(op, k): # tensor helper # ────────────────── single-pulse FULL ────────────────── -def run_simulation_single_pulse_full(freq, amp, tlist, init_freqs, interactions, - gamma, gamma_phi, pulse_ns, T_ramp_ns=1.0): +def run_simulation_single_pulse_full( + freq, + amp, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + pulse_ns, + T_ramp_ns=1.0, +): """ Returns pop(t), Sp(t), Sm(t) for one drive frequency & amplitude. freq, amp – as before (freq in “GHz” ≡ ns⁻¹ numerically). """ N = len(init_freqs) - sx1 = tensor(sigmax(), qeye(2)); sx2 = tensor(qeye(2), sigmax()) - sz1 = tensor(sigmaz(), qeye(2)); sz2 = tensor(qeye(2), sigmaz()) - sp1 = tensor(sigmap(), qeye(2)); sp2 = tensor(qeye(2), sigmap()) - sm1 = tensor(sigmam(), qeye(2)); sm2 = tensor(qeye(2), sigmam()) - + # sx1 = tensor(sigmax(), qeye(2)) + # sx2 = tensor(qeye(2), sigmax()) + # sz1 = tensor(sigmaz(), qeye(2)) + # sz2 = tensor(qeye(2), sigmaz()) + sp1 = tensor(sigmap(), qeye(2)) + sp2 = tensor(qeye(2), sigmap()) + # sm1 = tensor(sigmam(), qeye(2)) + # sm2 = tensor(qeye(2), sigmam()) sm, sx, _, sz = _tls_ops(N) - Sp = sum(s.dag() for s in sm) + # Sp = sum(s.dag() for s in sm) Sm = sum(sm) - Sx = sum(sx) + # Sx = sum(sx) - ## This was the original thing... + # This was the original thing... # T_ramp = 1.0 # def drive_coeff(t, args): # if t < T_ramp: @@ -70,12 +88,18 @@ def drive_coeff(t, args): if t <= pulse_ns: return amp * np.cos(freq * t) return 0.0 - - H_tls = sum(init_freqs[j]* sz[j] / 2 for j in range(N)) - # H_int = sum(interactions[i, j] * (sx[i] * sx[j]) for i in range(N_tls) for j in range(i)) # Dipole-Dipole Interaction - H_int = sum(interactions[i, j] for i in range(N) for j in range(i)) # Dipole-Dipole Interaction - + H_tls = sum(init_freqs[j] * sz[j] / 2 for j in range(N)) + + # Dipole-Dipole Interaction + # H_int = sum( + # interactions[i, j] * (sx[i] * sx[j]) + # for i in range(N_tls) for j in range(i) + # ) + H_int = sum( + interactions[i, j] for i in range(N) for j in range(i) + ) # Dipole-Dipole Interaction + H_drive = [[sum(sx), drive_coeff]] H = [H_tls + H_int] + H_drive @@ -86,7 +110,7 @@ def drive_coeff(t, args): if gamma_phi and gamma_phi > 0.0: # per-TLS pure dephasing (L = sqrt(gamma_phi/2) * sz_k) for k in range(N): - c_ops.append(np.sqrt(0.5*gamma_phi) * sz[k]) + c_ops.append(np.sqrt(0.5 * gamma_phi) * sz[k]) # ---------------------- Initial State ---------------------- H_static = H_tls + H_int @@ -94,12 +118,11 @@ def drive_coeff(t, args): psi0 = evecs[0] # Ground state # -------- Solve for and -------- - e_ops = [sp1, sp2, Sm.dag()*Sm] + e_ops = [sp1, sp2, Sm.dag() * Sm] # QuTiP-5 API compat: e_ops passed by keyword (now keyword-only); progress_bar # (a display flag only) moved into options. Physics unchanged. - result = mesolve(H, psi0, tlist, c_ops, e_ops=e_ops, - options={"progress_bar": None}) + result = mesolve(H, psi0, tlist, c_ops, e_ops=e_ops, options={"progress_bar": None}) Sp1_tot = result.expect[0] Sp2_tot = result.expect[1] @@ -107,53 +130,82 @@ def drive_coeff(t, args): return Sp1_tot, Sp2_tot, np.real(expec_pop) + # ────────────────── original thin wrapper ────────────── -def run_simulation_single_pulse(freq, amp, tlist, init_freqs, interactions, - gamma, gamma_phi, pulse_ns, T_ramp_ns=1.0): +def run_simulation_single_pulse( + freq, + amp, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + pulse_ns, + T_ramp_ns=1.0, +): sp1, sp2, pop = run_simulation_single_pulse_full( - freq, amp, tlist, init_freqs, interactions, - gamma, gamma_phi, pulse_ns, T_ramp_ns + freq, + amp, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + pulse_ns, + T_ramp_ns, ) return sp1, sp2, pop # ────────────────── two-pulse (unchanged) ────────────── -def run_simulation_double_pulse(freq, amp1, amp2, - pulse1_ns, gap_ns, pulse2_ns, - tlist, init_freqs, interactions, - gamma, gamma_phi, T_ramp_ns=1.0): +def run_simulation_double_pulse( + freq, + amp1, + amp2, + pulse1_ns, + gap_ns, + pulse2_ns, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + T_ramp_ns=1.0, +): sm, sx, _, sz = _tls_ops(len(init_freqs)) Sx, Sm = sum(sx), sum(sm) - H_tls = sum(init_freqs[j]*sz[j]/2 for j in range(len(init_freqs))) + H_tls = sum(init_freqs[j] * sz[j] / 2 for j in range(len(init_freqs))) t_start_p2 = pulse1_ns + gap_ns + def drive_coeff(t, _=None): if t < T_ramp_ns: - env = 0.5*(1-np.cos(np.pi*t/T_ramp_ns)) - return amp1*env*np.cos(freq*t) + env = 0.5 * (1 - np.cos(np.pi * t / T_ramp_ns)) + return amp1 * env * np.cos(freq * t) elif t < pulse1_ns: - return amp1*np.cos(freq*t) + return amp1 * np.cos(freq * t) if t < t_start_p2: return 0.0 t2 = t - t_start_p2 if t2 < T_ramp_ns: - env = 0.5*(1-np.cos(np.pi*t2/T_ramp_ns)) - return amp2*env*np.cos(freq*t) + env = 0.5 * (1 - np.cos(np.pi * t2 / T_ramp_ns)) + return amp2 * env * np.cos(freq * t) elif t2 < pulse2_ns: - return amp2*np.cos(freq*t) + return amp2 * np.cos(freq * t) return 0.0 H = [H_tls + interactions, [Sx, drive_coeff]] psi0 = (H_tls + interactions).groundstate()[1] - c_ops = [np.sqrt(gamma)*Sm] + c_ops = [np.sqrt(gamma) * Sm] if gamma_phi: _, _, _, szs = _tls_ops(len(init_freqs)) - c_ops += [np.sqrt(gamma_phi)*z for z in szs] + c_ops += [np.sqrt(gamma_phi) * z for z in szs] return np.real_if_close( # QuTiP-5 API compat: e_ops keyword, progress_bar into options. Physics unchanged. - mesolve(H, psi0, tlist, c_ops, e_ops=[Sm.dag()*Sm], - options={"progress_bar": None}).expect[0] + mesolve( + H, psi0, tlist, c_ops, e_ops=[Sm.dag() * Sm], options={"progress_bar": None} + ).expect[0] ) @@ -172,45 +224,54 @@ def drive_coeff(t, _=None): # ) # return H_int + # -------------------------------------------------- -def build_spin_spin_interactions_random_distribution(N_tls, J_min, J_max, - alpha_x=1.0, alpha_y=1.0, alpha_z=0.5): - +def build_spin_spin_interactions_random_distribution( + N_tls, J_min, J_max, alpha_x=1.0, alpha_y=1.0, alpha_z=0.5 +): + I2 = qeye(2) sm_single = destroy(2) sx_single = sm_single + sm_single.dag() - sy_single = -1j*(sm_single - sm_single.dag()) - sz_single = 2*sm_single.dag()*sm_single - I2 - + sy_single = -1j * (sm_single - sm_single.dag()) + sz_single = 2 * sm_single.dag() * sm_single - I2 + def tensor_op(op, k, N): - return tensor([I2]*k + [op] + [I2]*(N - k - 1)) - + return tensor([I2] * k + [op] + [I2] * (N - k - 1)) + sx_ops = [tensor_op(sx_single, i, N_tls) for i in range(N_tls)] sy_ops = [tensor_op(sy_single, i, N_tls) for i in range(N_tls)] sz_ops = [tensor_op(sz_single, i, N_tls) for i in range(N_tls)] - + J_mat = np.zeros((N_tls, N_tls), dtype=np.complex128) for i in range(N_tls): - for j in range(i+1, N_tls): + for j in range(i + 1, N_tls): np.random.seed(42) J_rand = np.random.uniform(J_min, J_max) - J_mat[i,j] = J_rand - J_mat[j,i] = J_rand - + J_mat[i, j] = J_rand + J_mat[j, i] = J_rand + H_int = 0 for i in range(N_tls): - for j in range(i+1, N_tls): - Jij = J_mat[i,j] - H_int += alpha_x*Jij*(sx_ops[i]*sx_ops[j]) - H_int += alpha_y*Jij*(sy_ops[i]*sy_ops[j]) - H_int += alpha_z*Jij*(sz_ops[i]*sz_ops[j]) + for j in range(i + 1, N_tls): + Jij = J_mat[i, j] + H_int += alpha_x * Jij * (sx_ops[i] * sx_ops[j]) + H_int += alpha_y * Jij * (sy_ops[i] * sy_ops[j]) + H_int += alpha_z * Jij * (sz_ops[i] * sz_ops[j]) return H_int # ───────── legacy wrapper (unchanged) ───────── -def run_simulation_for_frequency(freq, tlist, init_freqs, interactions, - gamma, gamma_phi, drive_ampl, pulse_duration): +def run_simulation_for_frequency( + freq, tlist, init_freqs, interactions, gamma, gamma_phi, drive_ampl, pulse_duration +): return run_simulation_single_pulse( - freq, drive_ampl, tlist, init_freqs, interactions, - gamma, gamma_phi, pulse_duration + freq, + drive_ampl, + tlist, + init_freqs, + interactions, + gamma, + gamma_phi, + pulse_duration, ) diff --git a/src/dataset_creation_scripts/11_temperature_sweep_warmup/experiment_11_dataset_creation.py b/src/dataset_creation_scripts/11_temperature_sweep_warmup/experiment_11_dataset_creation.py index 8f5f048..3d96358 100644 --- a/src/dataset_creation_scripts/11_temperature_sweep_warmup/experiment_11_dataset_creation.py +++ b/src/dataset_creation_scripts/11_temperature_sweep_warmup/experiment_11_dataset_creation.py @@ -144,7 +144,6 @@ import numpy as np from tqdm import tqdm - SAMPLES = [("FM_Shipley", 0)] SAMPLE_MATERIAL = { "FM_Shipley": "Sapphire substrate with spin-coated Shipley 1813 photoresist" @@ -173,7 +172,7 @@ PULSE_WIDTH_TICKS = 308 AWG_CLOCK_MHZ = 9830.4 # pulse generator's clock -> 1 tick ~= 0.1017 ns ADC_CLOCK_MHZ = 552.96 # digitiser's clock -> 1 sample ~= 1.8084 ns -N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down +N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down N_FREQUENCIES = len(EXPECTED_FREQ_MHZ) # 401 PULSE_WIDTH_NS = (PULSE_WIDTH_TICKS / AWG_CLOCK_MHZ) * 1000 # ~= 31.33 ns @@ -566,19 +565,26 @@ def build_dataset( n_reps = len(reps) n_total_rows = n_reps * N_FREQUENCIES * N_TIME_SAMPLES - print(f"\nFormat : long (one row per time sample)") + print("\nFormat : long (one row per time sample)") print(f"Sample : {SAMPLE_NAMES[0]} (sapphire + Shipley 1813)") print(f"Reps found : {n_reps} ({min(reps)}-{max(reps)})") - print(f"Frequencies : {N_FREQUENCIES} ({FREQ_START_MHZ}-{FREQ_STOP_MHZ} MHz, " - f"step {FREQ_STEP_MHZ})") + print( + f"Frequencies : {N_FREQUENCIES} ({FREQ_START_MHZ}-{FREQ_STOP_MHZ} MHz, " + f"step {FREQ_STEP_MHZ})" + ) print(f"Time samples : {N_TIME_SAMPLES} per measurement") - print(f"Total rows : {n_reps} x {N_FREQUENCIES} x {N_TIME_SAMPLES} " - f"= {n_total_rows:,}") + print( + f"Total rows : {n_reps} x {N_FREQUENCIES} x {N_TIME_SAMPLES} " + f"= {n_total_rows:,}" + ) print(f"Columns : {COLUMN_NAMES}") - print(f"Memory needed : ~{n_total_rows*len(COLUMN_NAMES)*4/1024**3:.2f} GB " - f"for the array alone") - print(f"Note : single-pulse experiment; T_mxc_K may legitimately " - f"be NaN\n") + print( + f"Memory needed : ~{n_total_rows * len(COLUMN_NAMES) * 4 / 1024**3:.2f} GB " + f"for the array alone" + ) + print( + "Note : single-pulse experiment; T_mxc_K may legitimately be NaN\n" + ) # ── Pass 1: find the run start, so elapsed_s is relative to it ─────────── run_start = None @@ -594,7 +600,7 @@ def build_dataset( # ── Pass 2: build ──────────────────────────────────────────────────────── data_all = np.empty((n_total_rows, len(COLUMN_NAMES)), dtype=np.float32) row_idx = 0 - sample_id = SAMPLES[0][1] + # sample_id = SAMPLES[0][1] n_files = 0 for rep in tqdm(reps, desc="Reps"): @@ -623,14 +629,18 @@ def build_dataset( T_mxc = z["T_mxc"] if IQ_matrix.shape != (2, N_FREQUENCIES, N_TIME_SAMPLES): - print(f" WARNING: {filename} has shape {IQ_matrix.shape}, expected " - f"(2, {N_FREQUENCIES}, {N_TIME_SAMPLES}) — skipping") + print( + f" WARNING: {filename} has shape {IQ_matrix.shape}, expected " + f"(2, {N_FREQUENCIES}, {N_TIME_SAMPLES}) — skipping" + ) continue # The frequency axis is in the file; verify rather than assume. if not np.array_equal(freq_array, EXPECTED_FREQ_MHZ): - print(f" WARNING: {filename} pulse_freq_array does not match " - f"np.arange(3000, 5001, 5) — using the file's own axis") + print( + f" WARNING: {filename} pulse_freq_array does not match " + f"np.arange(3000, 5001, 5) — using the file's own axis" + ) elapsed = (time_stamps - run_start).astype("timedelta64[s]").astype(np.float64) @@ -677,28 +687,30 @@ def build_dataset( with open(cols_path, "w", encoding="utf-8") as fh: json.dump(build_column_info(reps), fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)})") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") # ── Sanity checks ──────────────────────────────────────────────────────── - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] if row_idx != n_total_rows: errors.append(f" [FAIL] Row count: got {row_idx:,}, expected {n_total_rows:,}") else: - print(f" [PASS] Row count: {row_idx:,} = {n_reps} reps x {N_FREQUENCIES} " - f"frequencies x {N_TIME_SAMPLES} samples") + print( + f" [PASS] Row count: {row_idx:,} = {n_reps} reps x {N_FREQUENCIES} " + f"frequencies x {N_TIME_SAMPLES} samples" + ) # I and Q must NEVER be NaN, even though the thermometry may be. n_nan_iq = int(np.isnan(data_final[:, [COL["I"], COL["Q"]]]).sum()) if n_nan_iq: errors.append(f" [FAIL] I/Q contain {n_nan_iq} NaN values") else: - print(f" [PASS] No NaN in I/Q") + print(" [PASS] No NaN in I/Q") # NaN in T_mxc_K is expected: the thermometer drops out as the fridge warms. # Report where, rather than treating it as an error. @@ -706,15 +718,23 @@ def build_dataset( rep_col = data_final[:, COL["rep_index"]] n_nan_t = int(np.isnan(t_mxc).sum()) nan_reps = np.unique(rep_col[np.isnan(t_mxc)]).astype(int).tolist() - print(f" [INFO] T_mxc_K NaN: {n_nan_t:,} rows " - f"({n_nan_t/len(data_final)*100:.1f}%), in reps " - f"{min(nan_reps) if nan_reps else '-'}-" - f"{max(nan_reps) if nan_reps else '-'}") + print( + f" [INFO] T_mxc_K NaN: {n_nan_t:,} rows " + f"({n_nan_t / len(data_final) * 100:.1f}%), in reps " + f"{min(nan_reps) if nan_reps else '-'}-" + f"{max(nan_reps) if nan_reps else '-'}" + ) - I_min, I_max = float(data_final[:, COL["I"]].min()), float(data_final[:, COL["I"]].max()) - Q_min, Q_max = float(data_final[:, COL["Q"]].min()), float(data_final[:, COL["Q"]].max()) - print(f" [INFO] I range: [{I_min:.1f}, {I_max:.1f}] " - f"Q range: [{Q_min:.1f}, {Q_max:.1f}]") + I_min, I_max = float(data_final[:, COL["I"]].min()), float( + data_final[:, COL["I"]].max() + ) + Q_min, Q_max = float(data_final[:, COL["Q"]].min()), float( + data_final[:, COL["Q"]].max() + ) + print( + f" [INFO] I range: [{I_min:.1f}, {I_max:.1f}] " + f"Q range: [{Q_min:.1f}, {Q_max:.1f}]" + ) for col_idx, name, expected in [ (COL["rep_index"], "reps", n_reps), @@ -732,12 +752,16 @@ def build_dataset( if n_at == 0: errors.append(" [FAIL] frequency_MHz == 3400 matches no rows after upcast") else: - print(f" [PASS] frequency_MHz exact under float64 upcast " - f"({n_at:,} rows at 3400 MHz)") + print( + f" [PASS] frequency_MHz exact under float64 upcast " + f"({n_at:,} rows at 3400 MHz)" + ) el = data_final[:, COL["elapsed_s"]] - print(f" [INFO] elapsed_s range: [{el.min():.0f}, {el.max():.0f}] s " - f"({(el.max()-el.min())/3600:.1f} h warm-up)") + print( + f" [INFO] elapsed_s range: [{el.min():.0f}, {el.max():.0f}] s " + f"({(el.max() - el.min()) / 3600:.1f} h warm-up)" + ) tm = t_mxc[~np.isnan(t_mxc)] print(f" [INFO] T_mxc_K range: [{tm.min():.4f}, {tm.max():.2f}] K") @@ -748,24 +772,27 @@ def build_dataset( if not same: errors.append(" [FAIL] Pickle round-trip: reloaded data differs") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly " - f"(NaN-aware)") + print(" [PASS] Pickle round-trip: reloaded array matches exactly (NaN-aware)") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") - print(f" (expect roughly {n_total_rows*85/1024**3:.1f} GB and several " - f"minutes)") + print( + f" (expect roughly {n_total_rows * 85 / 1024**3:.1f} GB and several " + f"minutes)" + ) write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) if __name__ == "__main__": diff --git a/src/dataset_creation_scripts/12_phase_control_temperature/experiment_12_dataset_creation.py b/src/dataset_creation_scripts/12_phase_control_temperature/experiment_12_dataset_creation.py index 0440760..7eb1a78 100644 --- a/src/dataset_creation_scripts/12_phase_control_temperature/experiment_12_dataset_creation.py +++ b/src/dataset_creation_scripts/12_phase_control_temperature/experiment_12_dataset_creation.py @@ -171,7 +171,7 @@ PULSE_WIDTH_TICKS = 308 AWG_CLOCK_MHZ = 9830.4 # pulse generator's clock -> 1 tick ~= 0.1017 ns ADC_CLOCK_MHZ = 552.96 # digitiser's clock -> 1 sample ~= 1.8084 ns -N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down +N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down PULSE_WIDTH_NS = (PULSE_WIDTH_TICKS / AWG_CLOCK_MHZ) * 1000 # ~= 31.33 ns # The clock along each recorded trace. The digitiser samples every ~1.8 ns, so @@ -315,13 +315,11 @@ def build_metadata(campaign_start_iso: str = None, run_starts: dict = None) -> d "n": len(FREQUENCY_LIST_MHZ), "unit": "MHz", "values": FREQUENCY_LIST_MHZ, - }, "inter_pulse_spacing": { "n": len(SPACING_LIST_NS), "unit": "ns", "values": SPACING_LIST_NS, - }, "pulse_phase": { "n": len(EXPECTED_PHASE_DEG), @@ -556,15 +554,17 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): n_phases = len(EXPECTED_PHASE_DEG) n_total_rows = n_temps * n_freqs * n_spacings * n_phases * N_TIME_SAMPLES - print(f"\nFormat : long (one row per time sample)") + print("\nFormat : long (one row per time sample)") print(f"Sample : {SAMPLE_NAMES[0]} (sapphire + Shipley 1813)") print(f"Temperatures : {n_temps} ({[t for t, _ in TEMPERATURE_RUNS]} mK)") print(f"Frequencies : {n_freqs} ({FREQUENCY_LIST_MHZ} MHz)") print(f"Spacings : {n_spacings} ({SPACING_LIST_NS} ns)") print(f"Phases : {n_phases} (0-360 deg, step 3)") print(f"Time samples : {N_TIME_SAMPLES} per measurement") - print(f"Total rows : {n_temps} x {n_freqs} x {n_spacings} x {n_phases} x " - f"{N_TIME_SAMPLES} = {n_total_rows:,}") + print( + f"Total rows : {n_temps} x {n_freqs} x {n_spacings} x {n_phases} x " + f"{N_TIME_SAMPLES} = {n_total_rows:,}" + ) print(f"Columns : {COLUMN_NAMES}\n") # ── Pass 1: find the campaign start and each run's start ───────────────── @@ -597,7 +597,7 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): # ── Pass 2: build ──────────────────────────────────────────────────────── data_all = np.empty((n_total_rows, len(COLUMN_NAMES)), dtype=np.float32) row_idx = 0 - sample_id = SAMPLES[0][1] + # sample_id = SAMPLES[0][1] n_files = 0 for temperature_mk, exp_id in tqdm(TEMPERATURE_RUNS, desc="Temperatures"): @@ -620,13 +620,17 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): time_stamps = z["time_stamp_list"] if IQ_matrix.shape != (2, n_phases, N_TIME_SAMPLES): - print(f" WARNING: {filename} has shape {IQ_matrix.shape}, " - f"expected (2, {n_phases}, {N_TIME_SAMPLES}) — skipping") + print( + f" WARNING: {filename} has shape {IQ_matrix.shape}, " + f"expected (2, {n_phases}, {N_TIME_SAMPLES}) — skipping" + ) continue if not np.array_equal(phase_array, EXPECTED_PHASE_DEG): - print(f" WARNING: {filename} pulse_phase_array does not match " - f"np.arange(0, 361, 3) — using the file's own axis") + print( + f" WARNING: {filename} pulse_phase_array does not match " + f"np.arange(0, 361, 3) — using the file's own axis" + ) # Wall-clock time, i.e. when this reading was actually taken in # the lab. Stored as SECONDS SINCE THE CAMPAIGN STARTED (the 8 mK @@ -658,8 +662,12 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): data_all[row_idx:end_idx, COL["phase_deg"]] = phase_deg data_all[row_idx:end_idx, COL["timestamp_ns"]] = TIMESTAMPS_NS data_all[row_idx:end_idx, COL["elapsed_s"]] = elapsed[phase_idx] - data_all[row_idx:end_idx, COL["I"]] = I_matrix[phase_idx] # all 1000 - data_all[row_idx:end_idx, COL["Q"]] = Q_matrix[phase_idx] # all 1000 + data_all[row_idx:end_idx, COL["I"]] = I_matrix[ + phase_idx + ] # all 1000 + data_all[row_idx:end_idx, COL["Q"]] = Q_matrix[ + phase_idx + ] # all 1000 row_idx += N_TIME_SAMPLES @@ -685,33 +693,41 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): with open(cols_path, "w", encoding="utf-8") as fh: json.dump(build_column_info(), fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)})") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") # ── Sanity checks ──────────────────────────────────────────────────────── - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] if row_idx != n_total_rows: errors.append(f" [FAIL] Row count: got {row_idx:,}, expected {n_total_rows:,}") else: - print(f" [PASS] Row count: {row_idx:,} = {n_temps} temps x {n_freqs} freqs " - f"x {n_spacings} spacings x {n_phases} phases x {N_TIME_SAMPLES}") + print( + f" [PASS] Row count: {row_idx:,} = {n_temps} temps x {n_freqs} freqs " + f"x {n_spacings} spacings x {n_phases} phases x {N_TIME_SAMPLES}" + ) iq = data_final[:, [COL["I"], COL["Q"]]] n_nan = int(np.isnan(iq).sum()) if n_nan: errors.append(f" [FAIL] I/Q contain {n_nan} NaN values") else: - print(f" [PASS] No NaN in I/Q") + print(" [PASS] No NaN in I/Q") - I_min, I_max = float(data_final[:, COL["I"]].min()), float(data_final[:, COL["I"]].max()) - Q_min, Q_max = float(data_final[:, COL["Q"]].min()), float(data_final[:, COL["Q"]].max()) - print(f" [INFO] I range: [{I_min:.1f}, {I_max:.1f}] " - f"Q range: [{Q_min:.1f}, {Q_max:.1f}]") + I_min, I_max = float(data_final[:, COL["I"]].min()), float( + data_final[:, COL["I"]].max() + ) + Q_min, Q_max = float(data_final[:, COL["Q"]].min()), float( + data_final[:, COL["Q"]].max() + ) + print( + f" [INFO] I range: [{I_min:.1f}, {I_max:.1f}] " + f"Q range: [{Q_min:.1f}, {Q_max:.1f}]" + ) for col_idx, name, expected in [ (COL["temperature_mK"], "temperatures", n_temps), @@ -730,32 +746,38 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): if n_at_100 == 0: errors.append(" [FAIL] spacing_ns == 100 matches no rows after upcast") else: - print(f" [PASS] spacing_ns exact under float64 upcast " - f"({n_at_100:,} rows at 100 ns)") + print( + f" [PASS] spacing_ns exact under float64 upcast " + f"({n_at_100:,} rows at 100 ns)" + ) el = data_final[:, COL["elapsed_s"]] - print(f" [INFO] elapsed_s range: [{el.min():.0f}, {el.max():.0f}] s " - f"({(el.max()-el.min())/3600:.1f} h campaign)") + print( + f" [INFO] elapsed_s range: [{el.min():.0f}, {el.max():.0f}] s " + f"({(el.max() - el.min()) / 3600:.1f} h campaign)" + ) reloaded = load_pickle(pkl_path) if not np.array_equal(reloaded["data"], data_final): errors.append(" [FAIL] Pickle round-trip: reloaded data differs") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly") + print(" [PASS] Pickle round-trip: reloaded array matches exactly") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) if __name__ == "__main__": diff --git a/src/dataset_creation_scripts/13_thermal_cycle/experiment_13_dataset_creation.py b/src/dataset_creation_scripts/13_thermal_cycle/experiment_13_dataset_creation.py index c0026d4..4de6666 100644 --- a/src/dataset_creation_scripts/13_thermal_cycle/experiment_13_dataset_creation.py +++ b/src/dataset_creation_scripts/13_thermal_cycle/experiment_13_dataset_creation.py @@ -183,7 +183,7 @@ FREQUENCY_LIST_MHZ = np.arange(FREQ_START_MHZ, FREQ_STOP_MHZ, FREQ_STEP_MHZ) # 3000 ADC_CLOCK_MHZ = 552.96 # digitiser's clock -> 1 sample ~= 1.8084 ns -N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down +N_TIME_SAMPLES = 1000 # how many instants are recorded per ring-down N_FREQUENCIES = len(FREQUENCY_LIST_MHZ) # 3000 # The clock along each recorded trace. The digitiser samples every ~1.8 ns, so @@ -443,21 +443,25 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): n_cooldowns = len(COOLDOWNS) n_total_rows = n_cooldowns * N_FREQUENCIES * N_TIME_SAMPLES # 6,000,000 - print(f"\nFormat : long (one row per time sample)") + print("\nFormat : long (one row per time sample)") print(f"Sample : {SAMPLE_NAMES[0]} (sapphire + Shipley 1813)") print(f"Cooldowns : {n_cooldowns} (1 = before, 2 = after thermal cycle)") print(f"Temperature : ~{TEMPERATURE_MK} mK (nominal; not in the raw files)") - print(f"Frequencies : {N_FREQUENCIES} ({FREQUENCY_LIST_MHZ[0]}-" - f"{FREQUENCY_LIST_MHZ[-1]} MHz, step {FREQ_STEP_MHZ})") + print( + f"Frequencies : {N_FREQUENCIES} ({FREQUENCY_LIST_MHZ[0]}-" + f"{FREQUENCY_LIST_MHZ[-1]} MHz, step {FREQ_STEP_MHZ})" + ) print(f"Time samples : {N_TIME_SAMPLES} per measurement") - print(f"Total rows : {n_cooldowns} x {N_FREQUENCIES} x {N_TIME_SAMPLES} " - f"= {n_total_rows:,}") + print( + f"Total rows : {n_cooldowns} x {N_FREQUENCIES} x {N_TIME_SAMPLES} " + f"= {n_total_rows:,}" + ) print(f"Columns : {COLUMN_NAMES}") - print(f"Note : single-pulse experiment; pulse width is not recorded\n") + print("Note : single-pulse experiment; pulse width is not recorded\n") data_all = np.empty((n_total_rows, len(COLUMN_NAMES)), dtype=np.float32) row_idx = 0 - sample_id = SAMPLES[0][1] + # sample_id = SAMPLES[0][1] for cooldown_index, prefix, _desc in tqdm(COOLDOWNS, desc="Cooldowns"): @@ -471,8 +475,10 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): IQ_matrix = np.load(filepath) # expect (2, 3000, 1000) if IQ_matrix.shape != (2, N_FREQUENCIES, N_TIME_SAMPLES): - print(f" WARNING: {filename} has shape {IQ_matrix.shape}, expected " - f"(2, {N_FREQUENCIES}, {N_TIME_SAMPLES}) — skipping") + print( + f" WARNING: {filename} has shape {IQ_matrix.shape}, expected " + f"(2, {N_FREQUENCIES}, {N_TIME_SAMPLES}) — skipping" + ) continue I_matrix = IQ_matrix[0] # (3000, 1000) @@ -510,33 +516,41 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): with open(cols_path, "w", encoding="utf-8") as fh: json.dump(build_column_info(), fh, indent=2) - print(f"\nDone!") + print("\nDone!") print(f" Shape : ({row_idx:,}, {len(COLUMN_NAMES)})") - print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path)/1024**2:.1f} MB)") + print(f" Pickle : {pkl_path} ({os.path.getsize(pkl_path) / 1024**2:.1f} MB)") print(f" metadata : {meta_path}") print(f" column_info : {cols_path}") # ── Sanity checks ──────────────────────────────────────────────────────── - print(f"\nRunning sanity checks ...") + print("\nRunning sanity checks ...") errors = [] if row_idx != n_total_rows: errors.append(f" [FAIL] Row count: got {row_idx:,}, expected {n_total_rows:,}") else: - print(f" [PASS] Row count: {row_idx:,} = {n_cooldowns} cooldowns x " - f"{N_FREQUENCIES} frequencies x {N_TIME_SAMPLES} samples") + print( + f" [PASS] Row count: {row_idx:,} = {n_cooldowns} cooldowns x " + f"{N_FREQUENCIES} frequencies x {N_TIME_SAMPLES} samples" + ) iq = data_final[:, [COL["I"], COL["Q"]]] n_nan = int(np.isnan(iq).sum()) if n_nan: errors.append(f" [FAIL] I/Q contain {n_nan} NaN values") else: - print(f" [PASS] No NaN in I/Q") + print(" [PASS] No NaN in I/Q") - I_min, I_max = float(data_final[:, COL["I"]].min()), float(data_final[:, COL["I"]].max()) - Q_min, Q_max = float(data_final[:, COL["Q"]].min()), float(data_final[:, COL["Q"]].max()) - print(f" [INFO] I range: [{I_min:.1f}, {I_max:.1f}] " - f"Q range: [{Q_min:.1f}, {Q_max:.1f}]") + I_min, I_max = float(data_final[:, COL["I"]].min()), float( + data_final[:, COL["I"]].max() + ) + Q_min, Q_max = float(data_final[:, COL["Q"]].min()), float( + data_final[:, COL["Q"]].max() + ) + print( + f" [INFO] I range: [{I_min:.1f}, {I_max:.1f}] " + f"Q range: [{Q_min:.1f}, {Q_max:.1f}]" + ) for col_idx, name, expected in [ (COL["cooldown_index"], "cooldowns", n_cooldowns), @@ -555,8 +569,10 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): if n_at == 0: errors.append(" [FAIL] frequency_MHz == 3657 matches no rows after upcast") else: - print(f" [PASS] frequency_MHz exact under float64 upcast " - f"({n_at:,} rows at 3657 MHz)") + print( + f" [PASS] frequency_MHz exact under float64 upcast " + f"({n_at:,} rows at 3657 MHz)" + ) fr = data_final[:, COL["frequency_MHz"]] print(f" [INFO] frequency range: [{fr.min():.0f}, {fr.max():.0f}] MHz") @@ -565,21 +581,23 @@ def build_dataset(data_dir: str, output_dir: str, write_csv_file: bool = True): if not np.array_equal(reloaded["data"], data_final): errors.append(" [FAIL] Pickle round-trip: reloaded data differs") else: - print(f" [PASS] Pickle round-trip: reloaded array matches exactly") + print(" [PASS] Pickle round-trip: reloaded array matches exactly") if errors: print(f"\n {len(errors)} sanity check(s) FAILED:") for e in errors: print(e) else: - print(f"\n All sanity checks passed.") + print("\n All sanity checks passed.") # ── CSV ────────────────────────────────────────────────────────────────── if write_csv_file: print(f"\nWriting CSV to {csv_path} ...") write_csv(data_final, COLUMN_NAMES, csv_path) - print(f" CSV : {csv_path} " - f"({os.path.getsize(csv_path)/1024**2:.1f} MB, {row_idx:,} rows)") + print( + f" CSV : {csv_path} " + f"({os.path.getsize(csv_path) / 1024**2:.1f} MB, {row_idx:,} rows)" + ) if __name__ == "__main__": diff --git a/website/build_site.py b/website/build_site.py index 9a5d492..53e379e 100644 --- a/website/build_site.py +++ b/website/build_site.py @@ -1,16 +1,20 @@ #!/usr/bin/env python3 """Build the QTLS dataset portal: reads every experiment README and embeds it so each experiment card opens its own documentation in a reader panel.""" + import re, json, os HERE = os.path.dirname(os.path.abspath(__file__)) BASE = os.path.join(HERE, "..", "src", "dataset_creation_scripts") -OUT = os.path.join(HERE, "index.html") +OUT = os.path.join(HERE, "index.html") -SIMS = {6, 7, 8, 9, 10, 15, 16} # simulation experiments; the rest are measured +SIMS = {6, 7, 8, 9, 10, 15, 16} # simulation experiments; the rest are measured -folders = sorted(d for d in os.listdir(BASE) - if re.match(r"^\d\d_", d) and os.path.isdir(os.path.join(BASE, d))) +folders = sorted( + d + for d in os.listdir(BASE) + if re.match(r"^\d\d_", d) and os.path.isdir(os.path.join(BASE, d)) +) experiments = [] for folder in folders: @@ -55,32 +59,37 @@ if len(desc) > 185: desc = desc[:185].rsplit(" ", 1)[0] + "…" - experiments.append({ - "num": f"{num:02d}", - "title": title, - "desc": desc, - "badge": "sim" if num in SIMS else "meas", - "md": md, - }) - -readmes = {e["num"]: {"title": e["title"], "badge": e["badge"], "md": e["md"]} - for e in experiments} + experiments.append( + { + "num": f"{num:02d}", + "title": title, + "desc": desc, + "badge": "sim" if num in SIMS else "meas", + "md": md, + } + ) + +readmes = { + e["num"]: {"title": e["title"], "badge": e["badge"], "md": e["md"]} + for e in experiments +} # experiment cards cards = [] for e in experiments: label = "Simulated" if e["badge"] == "sim" else "Measured" cards.append( - f'''''') + """ + ) cards_html = "\n".join(cards) n_meas = sum(1 for e in experiments if e["badge"] == "meas") -n_sim = sum(1 for e in experiments if e["badge"] == "sim") +n_sim = sum(1 for e in experiments if e["badge"] == "sim") TEMPLATE = r""" @@ -514,12 +523,13 @@ """ -html = (TEMPLATE - .replace("__CARDS__", cards_html) - .replace("__READMES__", json.dumps(readmes)) - .replace("__NEXP__", str(len(experiments))) - .replace("__NMEAS__", str(n_meas)) - .replace("__NSIM__", str(n_sim))) +html = ( + TEMPLATE.replace("__CARDS__", cards_html) + .replace("__READMES__", json.dumps(readmes)) + .replace("__NEXP__", str(len(experiments))) + .replace("__NMEAS__", str(n_meas)) + .replace("__NSIM__", str(n_sim)) +) open(OUT, "w", encoding="utf-8").write(html) print("wrote", OUT, "-", len(experiments), "experiments,", len(html), "bytes")