GGN approximation formula defined

An approximated version of the GGN is implemented to reduce the computational time enabling fast multi-band transmission simulations

Change-Id: I2951a878aa04b5eb4a33ba86d626a788c4cbb100

docs/json.rst

@@ -523,6 +523,10 @@ See ``delta_power_range_db`` for more explaination.
 |                                             |           | ``ggn_spectrally_separated`` (see eq. 21    |
 |                                             |           | from `arXiv:1710.02225                      |
 |                                             |           | <https://arxiv.org/abs/1710.02225>`_).      |
+|                                             |           | ``ggn_approx`` (see eq. 24-25               |
+|                                             |           | from `jlt:9741324                           |
+|                                             |           | <https://eeexplore.ieee.org/document/       |
+|                                             |           | 9741324>`_).                                |
 +---------------------------------------------+-----------+---------------------------------------------+
 | ``dispersion_tolerance``                    | (number)  | Optional. Pure number. Tuning parameter for |
 |                                             |           | ggn model solution. Default value is 1.     |

gnpy/core/science_utils.py

@@ -12,7 +12,7 @@
 
 from numpy import interp, pi, zeros, cos, array, append, ones, exp, arange, sqrt, trapz, arcsinh, clip, abs, sum, \
     concatenate, flip, outer, inner, transpose, max, format_float_scientific, diag, sort, unique, argsort, cumprod, \
-    polyfit
+    polyfit, log, reshape, swapaxes, full, nan
 from logging import getLogger
 from scipy.constants import k, h
 from scipy.interpolate import interp1d
@@ -276,6 +276,7 @@ class NliSolver:
     List of implemented methods:
     'gn_model_analytic': eq. 120 from arXiv:1209.0394
     'ggn_spectrally_separated': eq. 21 from arXiv: 1710.02225
+    'ggn_approx': eq. 24-25 jlt:9741324
     """
 
     SPM_WEIGHT = (16.0 / 27.0)
@@ -324,6 +325,28 @@ def compute_nli(spectral_info: SpectralInformation, srs: StimulatedRamanScatteri
             g_nli = sum(g_nli, 1)
             g_nli = interp(spectral_info.frequency, cut_frequency, g_nli)
             nli = spectral_info.baud_rate * g_nli  # Local white noise
+        elif 'ggn_approx' in sim_params.nli_params.method:
+            if sim_params.nli_params.computed_channels is not None:
+                cut_indices = array(sim_params.nli_params.computed_channels) - 1
+            elif sim_params.nli_params.computed_number_of_channels is not None:
+                nb_ch_computed = sim_params.nli_params.computed_number_of_channels
+                nb_ch = len(spectral_info.channel_number)
+                cut_indices = array([round(i * (nb_ch - 1) / (nb_ch_computed - 1)) for i in range(0, nb_ch_computed)])
+            else:
+                cut_indices = array(spectral_info.channel_number) - 1
+
+            eta = NliSolver._ggn_approx(cut_indices, spectral_info, fiber, srs)
+
+            # Interpolation over the channels not indicated as computed channels in simulation parameters
+            cut_power = outer(spectral_info.signal[cut_indices], ones(spectral_info.number_of_channels))
+            cut_frequency = spectral_info.frequency[cut_indices]
+            pump_power = outer(ones(cut_indices.size), spectral_info.signal)
+            cut_baud_rate = outer(spectral_info.baud_rate[cut_indices], ones(spectral_info.number_of_channels))
+
+            g_nli = eta * cut_power * pump_power ** 2 / cut_baud_rate
+            g_nli = sum(g_nli, 1)
+            g_nli = interp(spectral_info.frequency, cut_frequency, g_nli)
+            nli = spectral_info.baud_rate * g_nli  # Local white noise
         else:
             raise ValueError(f'Method {sim_params.nli_params.method} not implemented.')
 
@@ -526,6 +549,89 @@ def _frequency_offset_threshold(beta2, symbol_rate):
         freq_offset_th = ((k_ref * delta_f_ref) * rs_ref * beta2_ref) / (beta2 * symbol_rate)
         return freq_offset_th
 
+    @staticmethod
+    def _ggn_approx(cut_indices, spectral_info: SpectralInformation, fiber, srs, spm_weight=SPM_WEIGHT,
+                    xpm_weight=XPM_WEIGHT):
+        """Computes the nonlinear interference power evaluated at the fiber input.
+        The method uses eq. 24-25 of https://ieeexplore.ieee.org/document/9741324
+        """
+        # Spectral Features
+        nch = spectral_info.number_of_channels
+        frequency = spectral_info.frequency
+        baud_rate = spectral_info.baud_rate
+        slot_width = spectral_info.slot_width
+        roll_off = spectral_info.roll_off
+        df = spectral_info.df + diag(full(nch, nan))
+
+        # Physical fiber parameters
+        alpha = fiber.alpha(frequency)
+        beta2 = fiber.beta2(frequency)
+        gamma = outer(fiber.gamma(frequency[cut_indices]), ones(nch))
+
+        identity = diag(ones(nch))
+        weight = spm_weight * identity + xpm_weight * (ones([nch, nch]) - identity)
+        weight = weight[cut_indices, :]
+
+        dispersion_tolerance = sim_params.nli_params.dispersion_tolerance
+        phase_shift_tolerance = sim_params.nli_params.phase_shift_tolerance
+        max_slot_width = max(slot_width)
+        max_beta2 = max(abs(beta2))
+        delta_z = sim_params.raman_params.result_spatial_resolution
+
+        # Approximation psi
+        loss_profile = srs.loss_profile[:nch]
+        z = srs.z
+        psi = NliSolver._approx_psi(df=df, frequency=frequency, beta2=beta2, baud_rate=baud_rate,
+                                    loss_profile=loss_profile, z=z)
+
+        # GGN for SPM
+        for cut_index in cut_indices:
+            dn = 0
+            cut_frequency = frequency[cut_index]
+            cut_baud_rate = baud_rate[cut_index]
+            cut_roll_off = roll_off[cut_index]
+            cut_beta2 = beta2[cut_index]
+            cut_alpha = alpha[cut_index]
+            k_tol = dispersion_tolerance * abs(cut_alpha)
+            phi_tol = phase_shift_tolerance / delta_z
+            f_cut_resolution = min(k_tol, phi_tol) / abs(max_beta2) / (4 * pi ** 2 * (1 + dn) * max_slot_width)
+            f_pump_resolution = min(k_tol, phi_tol) / abs(max_beta2) / (4 * pi ** 2 * max_slot_width)
+            psi[cut_index, cut_index] = NliSolver._generalized_psi(cut_frequency, cut_frequency, cut_baud_rate,
+                                                                   cut_roll_off, cut_frequency, cut_baud_rate,
+                                                                   cut_roll_off, f_cut_resolution, f_pump_resolution,
+                                                                   srs, cut_alpha, cut_beta2, 0, cut_frequency)
+        psi = psi[cut_indices, :]
+        cut_baud_rate = outer(baud_rate[cut_indices], ones(nch))
+        pump_baud_rate = outer(ones(cut_indices.size), baud_rate)
+
+        eta_cut_central_frequency = \
+            gamma ** 2 * weight * psi / (cut_baud_rate * pump_baud_rate ** 2)
+        eta = cut_baud_rate * eta_cut_central_frequency  # Local white noise
+
+        return eta
+
+    @staticmethod
+    def _approx_psi(df, frequency, baud_rate, beta2, loss_profile, z):
+        """Computes the approximated psi function similarly to the one used in the GN model.
+        The method uses eq. 25 of https://ieeexplore.ieee.org/document/9741324"""
+        pump_baud_rate = outer(ones(frequency.size), baud_rate)
+        cut_beta = outer(beta2, ones(frequency.size))
+        pump_beta = outer(ones(frequency.size), beta2)
+        delta_z = abs(z[:-1] - z[1:])
+
+        loss_lin = log(loss_profile)
+        pump_alpha = (loss_lin[:, 1:] - loss_lin[:, :-1]) / delta_z
+        leff = abs((loss_profile[:, 1:] - loss_profile[:, :-1]) / sqrt(abs(pump_alpha))) * pump_alpha / abs(pump_alpha)
+        leff = reshape(outer(leff, ones(z.size - 1)), newshape=[leff.shape[0], leff.shape[1], leff.shape[1]])
+        leff2 = leff * swapaxes(leff, 2, 1)
+        leff2 = sum(leff2, axis=(1, 2))
+        z_int = outer(ones(frequency.size), leff2)
+
+        delta_beta = (cut_beta + pump_beta) / 2
+        psi = z_int * pump_baud_rate / (4 * pi * abs(delta_beta * df))
+        return psi
+
+
 
 def estimate_nf_model(type_variety, gain_min, gain_max, nf_min, nf_max):
     if nf_min < -10:

tests/test_science_utils.py

@@ -10,6 +10,7 @@
 from pandas import read_csv
 from numpy.testing import assert_allclose
 from numpy import array
+from copy import deepcopy
 import pytest
 
 from gnpy.core.info import create_input_spectral_information, create_arbitrary_spectral_information
@@ -132,3 +133,26 @@ def test_fiber_lumped_losses_srs(set_sim_params):
     stimulated_raman_scattering = RamanSolver.calculate_attenuation_profile(spectral_info_input, fiber)
     power_profile = stimulated_raman_scattering.power_profile
     assert_allclose(power_profile, expected_power_profile, rtol=1e-3)
+
+
+@pytest.mark.usefixtures('set_sim_params')
+def test_nli_solver():
+    """Test the accuracy of NLI solver."""
+    fiber = Fiber(**load_json(TEST_DIR / 'data' / 'test_science_utils_fiber_config.json'))
+    fiber.ref_pch_in_dbm = 0.0
+    # fix grid spectral information generation
+    spectral_info_input = create_input_spectral_information(f_min=191.3e12, f_max=196.1e12, roll_off=0.15,
+                                                            baud_rate=32e9, spacing=50e9, tx_osnr=40.0,
+                                                            tx_power=1e-3)
+
+    SimParams.set_params(load_json(TEST_DIR / 'data' / 'sim_params.json'))
+    sim_params = SimParams()
+
+    spectral_info_output = fiber(deepcopy(spectral_info_input))
+
+    sim_params.nli_params.method = 'ggn_approx'
+    sim_params.raman_params.result_spatial_resolution = 10e3
+    spectral_info_output_approx = fiber(deepcopy(spectral_info_input))
+
+    assert_allclose(spectral_info_output.nli, spectral_info_output_approx.nli, rtol=1e-1)
+