From 1feb3daea4f00a2f6caac90577c2303cd721d56f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 21 Jul 2026 07:45:23 +0000 Subject: [PATCH 1/4] Vectorize coupling integral evaluation paths --- package/src/openflash/meem_engine.py | 14 +-- package/src/openflash/multi_equations.py | 140 +++++++++++++++-------- 2 files changed, 99 insertions(+), 55 deletions(-) diff --git a/package/src/openflash/meem_engine.py b/package/src/openflash/meem_engine.py index 6391f61a..35c7bbdb 100644 --- a/package/src/openflash/meem_engine.py +++ b/package/src/openflash/meem_engine.py @@ -121,9 +121,9 @@ def build_problem_cache(self, problem: 'MEEMProblem') -> ProblemCache: # 1. Pre-compute I_nm I_nm_vals_precomputed = [np.zeros((NMK[bd], NMK[bd+1]), dtype=complex) for bd in range(boundary_count - 1)] for bd in range(boundary_count - 1): - for n in range(NMK[bd]): - for m in range(NMK[bd + 1]): - I_nm_vals_precomputed[bd][n, m] = I_nm(n, m, bd, d, h) + n_modes = np.arange(NMK[bd])[:, None] + m_modes = np.arange(NMK[bd + 1])[None, :] + I_nm_vals_precomputed[bd] = I_nm(n_modes, m_modes, bd, d, h) cache._set_I_nm_vals(I_nm_vals_precomputed) # 2. Pre-defined Partials @@ -134,11 +134,9 @@ def build_problem_cache(self, problem: 'MEEMProblem') -> ProblemCache: # 3. Dynamic I_mk closure def _calculate_I_mk_vals(m0, m_k_arr, N_k_arr): - vals = np.zeros((NMK[boundary_count - 1], NMK[boundary_count]), dtype=complex) - for m in range(NMK[boundary_count - 1]): - for k in range(NMK[boundary_count]): - vals[m, k] = I_mk(m, k, boundary_count - 1, d, m0, h, m_k_arr, N_k_arr) - return vals + m_modes = np.arange(NMK[boundary_count - 1])[:, None] + k_modes = np.arange(NMK[boundary_count])[None, :] + return I_mk(m_modes, k_modes, boundary_count - 1, d, m0, h, m_k_arr, N_k_arr) # 4. Integration constants int_R1_store, int_R2_store = {}, {} diff --git a/package/src/openflash/multi_equations.py b/package/src/openflash/multi_equations.py index e683b14e..27a0bf2f 100644 --- a/package/src/openflash/multi_equations.py +++ b/package/src/openflash/multi_equations.py @@ -72,60 +72,106 @@ def m_k(NMK, m0, h): ############################################# # vertical eigenvector coupling computation +def _prepare_mode_grids(row_mode, col_mode): + row_arr = np.asarray(row_mode) + col_arr = np.asarray(col_mode) + scalar_input = row_arr.ndim == 0 and col_arr.ndim == 0 + + if scalar_input: + return row_arr.reshape(1, 1), col_arr.reshape(1, 1), True + + if row_arr.ndim <= 1 and col_arr.ndim <= 1: + row_grid, col_grid = np.meshgrid(row_arr.reshape(-1), col_arr.reshape(-1), indexing="ij") + else: + row_grid, col_grid = np.broadcast_arrays(row_arr, col_arr) + + return row_grid, col_grid, False + + def I_nm(n, m, i, d, h): # coupling integral for two i-type regions dj = max(d[i], d[i+1]) # integration bounds at -h and -d - if n == 0 and m == 0: - return h - dj - lambda1 = lambda_ni(n, i, h, d) - lambda2 = lambda_ni(m, i + 1, h, d) - if n == 0 and m >= 1: - if dj == d[i+1]: - return 0 - else: - return sqrt(2) * sin(lambda2 * (h - dj)) / lambda2 - if n >= 1 and m == 0: - if dj == d[i]: - return 0 - else: - return sqrt(2) * sin(lambda1 * (h - dj)) / lambda1 - else: - frac1 = sin((lambda1 + lambda2)*(h-dj))/(lambda1 + lambda2) - if lambda1 == lambda2: - frac2 = (h - dj) - else: - frac2 = sin((lambda1 - lambda2)*(h-dj))/(lambda1 - lambda2) - return frac1 + frac2 + n_grid, m_grid, scalar_input = _prepare_mode_grids(n, m) + delta = h - dj + out = np.zeros(n_grid.shape, dtype=float) + + mask_00 = (n_grid == 0) & (m_grid == 0) + out[mask_00] = delta + + mask_n0_mpos = (n_grid == 0) & (m_grid >= 1) + if np.any(mask_n0_mpos) and dj != d[i+1]: + lambda2 = lambda_ni(m_grid[mask_n0_mpos], i + 1, h, d) + out[mask_n0_mpos] = sqrt(2) * sin(lambda2 * delta) / lambda2 + + mask_npos_m0 = (n_grid >= 1) & (m_grid == 0) + if np.any(mask_npos_m0) and dj != d[i]: + lambda1 = lambda_ni(n_grid[mask_npos_m0], i, h, d) + out[mask_npos_m0] = sqrt(2) * sin(lambda1 * delta) / lambda1 + + mask_general = ~(mask_00 | mask_n0_mpos | mask_npos_m0) + if np.any(mask_general): + lambda1 = lambda_ni(n_grid[mask_general], i, h, d) + lambda2 = lambda_ni(m_grid[mask_general], i + 1, h, d) + frac1 = sin((lambda1 + lambda2) * delta) / (lambda1 + lambda2) + frac2 = np.empty_like(frac1, dtype=float) + eq_mask = np.isclose(lambda1, lambda2) + frac2[eq_mask] = delta + neq_mask = ~eq_mask + if np.any(neq_mask): + frac2[neq_mask] = sin((lambda1[neq_mask] - lambda2[neq_mask]) * delta) / (lambda1[neq_mask] - lambda2[neq_mask]) + out[mask_general] = frac1 + frac2 + + return out.item() if scalar_input else out # REVISED I_mk to accept m_k_arr and N_k_arr def I_mk(m, k, i, d, m0, h, m_k_arr, N_k_arr): # coupling integral for i and e-type regions - # Use the pre-computed array - local_m_k_k = m_k_arr[k] # Access directly from array - + m_grid, k_grid, scalar_input = _prepare_mode_grids(m, k) + m_grid = m_grid.astype(int, copy=False) + k_grid = k_grid.astype(int, copy=False) dj = d[i] - if m == 0 and k == 0: - if m0 == inf: return 0 - elif m0 * h < M0_H_THRESH: - return (1/sqrt(N_k_arr[0])) * sinh(m0 * (h - dj)) / m0 # Use N_k_arr[0] - else: # high m0h approximation - return sqrt(2 * h / m0) * (exp(- m0 * dj) - exp(m0 * dj - 2 * m0 * h)) - if m == 0 and k >= 1: - return (1/sqrt(N_k_arr[k])) * sin(local_m_k_k * (h - dj)) / local_m_k_k # Use N_k_arr[k] - if m >= 1 and k == 0: - if m0 == inf: return 0 - elif m0 * h < M0_H_THRESH: - num = (-1)**m * sqrt(2) * (1/sqrt(N_k_arr[0])) * m0 * sinh(m0 * (h - dj)) # Use N_k_arr[0] - else: # high m0h approximation - num = (-1)**m * 2 * sqrt(h * m0 ** 3) *(exp(- m0 * dj) - exp(m0 * dj - 2 * m0 * h)) - denom = (m0**2 + lambda_ni(m, i, h, d) **2) - return num/denom - else: - lambda1 = lambda_ni(m, i, h, d) - if abs(local_m_k_k) == lambda1: - return sqrt(2/N_k_arr[k]) * (h - dj)/2 + delta = h - dj + out = np.zeros(m_grid.shape, dtype=float) + + mask_m0_k0 = (m_grid == 0) & (k_grid == 0) + if np.any(mask_m0_k0) and m0 != inf: + if m0 * h < M0_H_THRESH: + out[mask_m0_k0] = (1 / sqrt(N_k_arr[0])) * sinh(m0 * delta) / m0 + else: + out[mask_m0_k0] = sqrt(2 * h / m0) * (exp(-m0 * dj) - exp(m0 * dj - 2 * m0 * h)) + + mask_m0_kpos = (m_grid == 0) & (k_grid >= 1) + if np.any(mask_m0_kpos): + k_local = k_grid[mask_m0_kpos] + local_mk = m_k_arr[k_local] + out[mask_m0_kpos] = (1 / sqrt(N_k_arr[k_local])) * sin(local_mk * delta) / local_mk + + mask_mpos_k0 = (m_grid >= 1) & (k_grid == 0) + if np.any(mask_mpos_k0) and m0 != inf: + m_local = m_grid[mask_mpos_k0] + if m0 * h < M0_H_THRESH: + num = ((-1) ** m_local) * sqrt(2) * (1 / sqrt(N_k_arr[0])) * m0 * sinh(m0 * delta) else: - frac1 = sin((local_m_k_k + lambda1)*(h-dj))/(local_m_k_k + lambda1) - frac2 = sin((local_m_k_k - lambda1)*(h-dj))/(local_m_k_k - lambda1) - return sqrt(2/N_k_arr[k]) * (frac1 + frac2)/2 # Use N_k_arr[k] + num = ((-1) ** m_local) * 2 * sqrt(h * m0 ** 3) * (exp(-m0 * dj) - exp(m0 * dj - 2 * m0 * h)) + denom = m0**2 + lambda_ni(m_local, i, h, d) ** 2 + out[mask_mpos_k0] = num / denom + + mask_general = ~(mask_m0_k0 | mask_m0_kpos | mask_mpos_k0) + if np.any(mask_general): + m_local = m_grid[mask_general] + k_local = k_grid[mask_general] + lambda1 = lambda_ni(m_local, i, h, d) + local_mk = m_k_arr[k_local] + norm = sqrt(2 / N_k_arr[k_local]) / 2 + eq_mask = np.isclose(np.abs(local_mk), lambda1) + general_vals = np.empty_like(lambda1, dtype=float) + general_vals[eq_mask] = norm[eq_mask] * delta + neq_mask = ~eq_mask + if np.any(neq_mask): + frac1 = sin((local_mk[neq_mask] + lambda1[neq_mask]) * delta) / (local_mk[neq_mask] + lambda1[neq_mask]) + frac2 = sin((local_mk[neq_mask] - lambda1[neq_mask]) * delta) / (local_mk[neq_mask] - lambda1[neq_mask]) + general_vals[neq_mask] = norm[neq_mask] * (frac1 + frac2) + out[mask_general] = general_vals + + return out.item() if scalar_input else out ############################################# # b-vector computation From 574d659469c39690c3e657677760df8237e0dbed Mon Sep 17 00:00:00 2001 From: Becca <32020426+rebeccamccabe@users.noreply.github.com> Date: Mon, 27 Jul 2026 11:36:54 -0400 Subject: [PATCH 2/4] Add env var for PR title to prevent code marks being interpreted as command Set PR_TITLE environment variable for base version computation. --- .github/workflows/pr-tagging.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/pr-tagging.yml b/.github/workflows/pr-tagging.yml index 29dff58f..5bc0c9fe 100644 --- a/.github/workflows/pr-tagging.yml +++ b/.github/workflows/pr-tagging.yml @@ -33,8 +33,9 @@ jobs: - name: Compute base version (bump patch/minor/major from PR title) id: compute_base + env: + PR_TITLE: ${{ github.event.pull_request.title }} run: | - PR_TITLE="${{ github.event.pull_request.title || '' }}" base=$(bash scripts/compute_base.sh "$PR_TITLE") echo "base=$base" >> $GITHUB_OUTPUT From 920f560b63d7aab040846133c6cffccdc0616767 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 27 Jul 2026 16:29:08 +0000 Subject: [PATCH 3/4] Make zero assignments explicit in vectorized masks --- package/src/openflash/multi_equations.py | 45 +++++++++++++++--------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/package/src/openflash/multi_equations.py b/package/src/openflash/multi_equations.py index 87bd1cba..c1886e38 100644 --- a/package/src/openflash/multi_equations.py +++ b/package/src/openflash/multi_equations.py @@ -92,20 +92,26 @@ def I_nm(n, m, i, d, h): # coupling integral for two i-type regions dj = max(d[i], d[i+1]) # integration bounds at -h and -d n_grid, m_grid, scalar_input = _prepare_mode_grids(n, m) delta = h - dj - out = np.zeros(n_grid.shape, dtype=float) + out = np.full(n_grid.shape, np.nan, dtype=float) mask_00 = (n_grid == 0) & (m_grid == 0) out[mask_00] = delta mask_n0_mpos = (n_grid == 0) & (m_grid >= 1) - if np.any(mask_n0_mpos) and dj != d[i+1]: - lambda2 = lambda_ni(m_grid[mask_n0_mpos], i + 1, h, d) - out[mask_n0_mpos] = sqrt(2) * sin(lambda2 * delta) / lambda2 + if np.any(mask_n0_mpos): + if dj != d[i+1]: + lambda2 = lambda_ni(m_grid[mask_n0_mpos], i + 1, h, d) + out[mask_n0_mpos] = sqrt(2) * sin(lambda2 * delta) / lambda2 + else: + out[mask_n0_mpos] = 0 mask_npos_m0 = (n_grid >= 1) & (m_grid == 0) - if np.any(mask_npos_m0) and dj != d[i]: - lambda1 = lambda_ni(n_grid[mask_npos_m0], i, h, d) - out[mask_npos_m0] = sqrt(2) * sin(lambda1 * delta) / lambda1 + if np.any(mask_npos_m0): + if dj != d[i]: + lambda1 = lambda_ni(n_grid[mask_npos_m0], i, h, d) + out[mask_npos_m0] = sqrt(2) * sin(lambda1 * delta) / lambda1 + else: + out[mask_npos_m0] = 0 mask_general = ~(mask_00 | mask_n0_mpos | mask_npos_m0) if np.any(mask_general): @@ -129,11 +135,13 @@ def I_mk(m, k, i, d, m0, h, m_k_arr, N_k_arr): # coupling integral for i and e-t k_grid = k_grid.astype(int, copy=False) dj = d[i] delta = h - dj - out = np.zeros(m_grid.shape, dtype=float) + out = np.full(m_grid.shape, np.nan, dtype=float) mask_m0_k0 = (m_grid == 0) & (k_grid == 0) - if np.any(mask_m0_k0) and m0 != inf: - if m0 * h < M0_H_THRESH: + if np.any(mask_m0_k0): + if m0 == inf: + out[mask_m0_k0] = 0 + elif m0 * h < M0_H_THRESH: out[mask_m0_k0] = (1 / sqrt(N_k_arr[0])) * sinh(m0 * delta) / m0 else: out[mask_m0_k0] = sqrt(2 * h / m0) * (exp(-m0 * dj) - exp(m0 * dj - 2 * m0 * h)) @@ -145,14 +153,17 @@ def I_mk(m, k, i, d, m0, h, m_k_arr, N_k_arr): # coupling integral for i and e-t out[mask_m0_kpos] = (1 / sqrt(N_k_arr[k_local])) * sin(local_mk * delta) / local_mk mask_mpos_k0 = (m_grid >= 1) & (k_grid == 0) - if np.any(mask_mpos_k0) and m0 != inf: - m_local = m_grid[mask_mpos_k0] - if m0 * h < M0_H_THRESH: - num = ((-1) ** m_local) * sqrt(2) * (1 / sqrt(N_k_arr[0])) * m0 * sinh(m0 * delta) + if np.any(mask_mpos_k0): + if m0 == inf: + out[mask_mpos_k0] = 0 else: - num = ((-1) ** m_local) * 2 * sqrt(h * m0 ** 3) * (exp(-m0 * dj) - exp(m0 * dj - 2 * m0 * h)) - denom = m0**2 + lambda_ni(m_local, i, h, d) ** 2 - out[mask_mpos_k0] = num / denom + m_local = m_grid[mask_mpos_k0] + if m0 * h < M0_H_THRESH: + num = ((-1) ** m_local) * sqrt(2) * (1 / sqrt(N_k_arr[0])) * m0 * sinh(m0 * delta) + else: + num = ((-1) ** m_local) * 2 * sqrt(h * m0 ** 3) * (exp(-m0 * dj) - exp(m0 * dj - 2 * m0 * h)) + denom = m0**2 + lambda_ni(m_local, i, h, d) ** 2 + out[mask_mpos_k0] = num / denom mask_general = ~(mask_m0_k0 | mask_m0_kpos | mask_mpos_k0) if np.any(mask_general): From f36dbfb83d5247d62b7d0a09ba68a12de4b1a115 Mon Sep 17 00:00:00 2001 From: Becca <32020426+rebeccamccabe@users.noreply.github.com> Date: Mon, 27 Jul 2026 12:56:57 -0400 Subject: [PATCH 4/4] Restore lost comments about high m0h Co-authored-by: Becca <32020426+rebeccamccabe@users.noreply.github.com> --- package/src/openflash/multi_equations.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/src/openflash/multi_equations.py b/package/src/openflash/multi_equations.py index c1886e38..b5c1a01d 100644 --- a/package/src/openflash/multi_equations.py +++ b/package/src/openflash/multi_equations.py @@ -143,7 +143,7 @@ def I_mk(m, k, i, d, m0, h, m_k_arr, N_k_arr): # coupling integral for i and e-t out[mask_m0_k0] = 0 elif m0 * h < M0_H_THRESH: out[mask_m0_k0] = (1 / sqrt(N_k_arr[0])) * sinh(m0 * delta) / m0 - else: + else: # high m0h approximation out[mask_m0_k0] = sqrt(2 * h / m0) * (exp(-m0 * dj) - exp(m0 * dj - 2 * m0 * h)) mask_m0_kpos = (m_grid == 0) & (k_grid >= 1) @@ -160,7 +160,7 @@ def I_mk(m, k, i, d, m0, h, m_k_arr, N_k_arr): # coupling integral for i and e-t m_local = m_grid[mask_mpos_k0] if m0 * h < M0_H_THRESH: num = ((-1) ** m_local) * sqrt(2) * (1 / sqrt(N_k_arr[0])) * m0 * sinh(m0 * delta) - else: + else: # high m0h approximation num = ((-1) ** m_local) * 2 * sqrt(h * m0 ** 3) * (exp(-m0 * dj) - exp(m0 * dj - 2 * m0 * h)) denom = m0**2 + lambda_ni(m_local, i, h, d) ** 2 out[mask_mpos_k0] = num / denom