Merge pull request #302 from Telecominfraproject/develop

Merge develop into master in preparation for the 1.8 release

Highlights for the upcoming release:

- Raman simulation
- automatic building of Docker images
- bugfixes, refactoring, CI and docs improvements

.docker-entry.sh

@@ -0,0 +1,3 @@
+#!/bin/bash
+cp -nr /oopt-gnpy/examples /shared
+exec "$@"

.docker-travis.sh

@@ -0,0 +1,47 @@
+#!/bin/bash
+
+set -e
+
+IMAGE_NAME=telecominfraproject/oopt-gnpy
+IMAGE_TAG=$(git describe --tags)
+
+ALREADY_FOUND=0
+docker pull ${IMAGE_NAME}:${IMAGE_TAG} && ALREADY_FOUND=1
+
+if [[ $ALREADY_FOUND == 0 ]]; then
+  docker build . -t ${IMAGE_NAME}
+  docker tag ${IMAGE_NAME} ${IMAGE_NAME}:${IMAGE_TAG}
+
+  # shared directory setup: do not clobber the real data
+  mkdir trash
+  cd trash
+  docker run -it --rm --volume $(pwd):/shared ${IMAGE_NAME} ./transmission_main_example.py
+else
+  echo "Image ${IMAGE_NAME}:${IMAGE_TAG} already available, will just update the other tags"
+fi
+
+docker images
+
+do_docker_login() {
+  echo "${DOCKER_PASSWORD}" | docker login -u "${DOCKER_USERNAME}" --password-stdin
+}
+
+if [[ "${TRAVIS_PULL_REQUEST}" == "false" ]]; then
+  if [[ "${TRAVIS_BRANCH}" == "develop" || "${TRAVIS_BRANCH}" == "docker" ]]; then
+    echo "Publishing latest"
+    docker tag ${IMAGE_NAME}:${IMAGE_TAG} ${IMAGE_NAME}:latest
+    do_docker_login
+    if [[ $ALREADY_FOUND == 0 ]]; then
+      docker push ${IMAGE_NAME}:${IMAGE_TAG}
+    fi
+    docker push ${IMAGE_NAME}:latest
+  elif [[ "${TRAVIS_BRANCH}" == "master" ]]; then
+    echo "Publishing stable"
+    docker tag ${IMAGE_NAME}:${IMAGE_TAG} ${IMAGE_NAME}:stable
+    do_docker_login
+    if [[ $ALREADY_FOUND == 0 ]]; then
+      docker push ${IMAGE_NAME}:${IMAGE_TAG}
+    fi
+    docker push ${IMAGE_NAME}:stable
+  fi
+fi

.travis.yml

@@ -1,15 +1,27 @@
 dist: xenial
 sudo: false
 language: python
+services: docker
 python:
   - "3.6"
   - "3.7"
-# command to install dependencies
+install: skip
-install:
+script:
   - python setup.py install
   - pip install pytest-cov rstcheck
-script:
   - pytest --cov-report=xml --cov=gnpy
   - rstcheck --ignore-roles cite --ignore-directives automodule --recursive --ignore-messages '(Duplicate explicit target name.*)' .
+  - ./examples/transmission_main_example.py
+  - ./examples/path_requests_run.py
+  - ./examples/transmission_main_example.py examples/raman_edfa_example_network.json --sim examples/sim_params.json --show-channels
+  - sphinx-build docs/ x-throwaway-location
 after_success:
   - bash <(curl -s https://codecov.io/bash)
+jobs:
+  include:
+    - stage: test
+      name: Docker image
+      script:
+        - git fetch --unshallow
+        - ./.docker-travis.sh
+        - docker images

Dockerfile

@@ -0,0 +1,7 @@
+FROM python:3.7-slim
+COPY . /oopt-gnpy
+WORKDIR /oopt-gnpy
+RUN python setup.py install
+WORKDIR /shared/examples
+ENTRYPOINT ["/oopt-gnpy/.docker-entry.sh"]
+CMD ["/bin/bash"]

Excel_userguide.rst

@@ -19,8 +19,8 @@ In order to work the excel file MUST contain at least 2 sheets:
 Nodes sheet
 -----------
 
-Nodes sheet contains seven columns.
+Nodes sheet contains nine columns.
-Each line represents a 'node' (ROADM site or an in line amplifier site ILA)::
+Each line represents a 'node' (ROADM site or an in line amplifier site ILA or a Fused)::
 
   City (Mandatory) ; State ; Country ; Region ; Latitude ; Longitude ; Type
 
@@ -38,6 +38,9 @@ Each line represents a 'node' (ROADM site or an in line amplifier site ILA)::
 
 - *Longitude*, *Latitude* are not mandatory. If filled they should contain numbers.
 
+- **Booster_restriction** and **Preamp_restriction** are not mandatory.
+  If used, they must contain one or several amplifier type_variety names separated by ' | '. This information is used to restrict types of amplifiers used in a ROADM node during autodesign. If a ROADM booster or preamp is already specified in the Eqpt sheet , the field is ignored. The field is also ignored if the node is not a ROADM node.
+
 **There MUST NOT be empty line(s) between two nodes lines**
 
 
@@ -166,6 +169,7 @@ This generates a text file meshTopologyExampleV2_eqt_sheet.txt  whose content ca
 - **amp type** is not mandatory. 
   If filled it must contain types listed in `eqpt_config.json <examples/eqpt_config.json>`_ in "Edfa" list "type_variety".
   If not filled it takes "std_medium_gain" as default value.
+  If filled with fused, a fused element with 0.0 dB loss will be placed instead of an amplifier. This might be used to avoid booster amplifier on a ROADM direction.
 
 - **amp_gain** is not mandatory. It is the value to be set on the amplifier (in dB).
   If not filled, it will be determined with design rules in the convert.py file.

README.rst

@@ -41,6 +41,33 @@ Branches and Tagged Releases
 How to Install
 --------------
 
+Using prebuilt Docker images
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Our `Docker images <https://hub.docker.com/r/telecominfraproject/oopt-gnpy>`_ contain everything needed to run all examples from this guide.
+Docker transparently fetches the image over the network upon first use.
+On Linux and Mac, run:
+
+
+.. code-block:: shell-session
+
+    $ docker run -it --rm --volume $(pwd):/shared telecominfraproject/oopt-gnpy
+    root@bea050f186f7:/shared/examples#
+
+On Windows, launch from Powershell as:
+
+.. code-block:: powershell
+
+    PS C:\> docker run -it --rm --volume ${PWD}:/shared telecominfraproject/oopt-gnpy
+    root@89784e577d44:/shared/examples#
+
+In both cases, a directory named ``examples/`` will appear in your current working directory.
+GNPy automaticallly populates it with example files from the current release.
+Remove that directory if you want to start from scratch.
+
+Using Python on your computer
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
    **Note**: `gnpy` supports Python 3 only. Python 2 is not supported.
    `gnpy` requires Python ≥3.6
 
@@ -105,7 +132,6 @@ executes without a ``ModuleNotFoundError``, you have successfully installed
 
     $ python -c 'import gnpy' # attempt to import gnpy
 
-    $ cd oopt-gnpy
     $ pytest                  # run tests
 
 Instructions for First Use
@@ -118,7 +144,7 @@ fully-functional programs.
 
     **Note**: *If you are a network operator or involved in route planning and
     optimization for your organization, please contact project maintainer Jan
-    Kundrát <jan.kundrat@telecominfraproject>. gnpy is looking for users with
+    Kundrát <jan.kundrat@telecominfraproject.com>. gnpy is looking for users with
     specific, delineated use cases to drive requirements for future
     development.*
 
@@ -417,10 +443,15 @@ existing parameters:
 +--------------------------+-----------+---------------------------------------------+
 | ``add_drop_osnr``        | (number)  | OSNR contribution from the add/drop ports   |
 +--------------------------+-----------+---------------------------------------------+
-| ``restrictions``         | (strings) | Authorized type_variety of amplifier for    |
+| ``restrictions``         | (dict of  | If non-empty, keys ``preamp_variety_list``  |
-|                          |           | booster or preamp.                          |
+|                          |  strings) | and ``booster_variety_list`` represent      |
-|                          |           | Listed type_variety MUST be defined in the  |
+|                          |           | list of ``type_variety`` amplifiers which   |
-|                          |           | Edfa catalog.                               |
+|                          |           | are allowed for auto-design within ROADM's  |
+|                          |           | line degrees.                               |
+|                          |           |                                             |
+|                          |           | If no booster should be placed on a degree, |
+|                          |           | insert a ``Fused`` node on the degree       |
+|                          |           | output.                                     |
 +--------------------------+-----------+---------------------------------------------+
 
 The ``SpectralInformation`` object can be configured as follows. The user can
@@ -470,6 +501,17 @@ Launch power can be overridden by using the ``--power`` argument.
 Spectrum information is not yet parametrized but can be modified directly in the ``eqpt_config.json`` (via the ``SpectralInformation`` -SI- structure) to accommodate any baud rate or spacing.
 The number of channel is computed based on ``spacing`` and ``f_min``, ``f_max`` values.
 
+An experimental support for Raman amplification is available:
+
+.. code-block:: shell
+
+     $ ./examples/transmission_main_example.py \
+       examples/raman_edfa_example_network.json \
+       --sim examples/sim_params.json --show-channels
+
+Configuration of Raman pumps (their frequencies, power and pumping direction) is done via the `RamanFiber element in the network topology <examples/raman_edfa_example_network.json>`_.
+General numeric parameters for simulaiton control are provided in the `examples/sim_params.json <examples/sim_params.json>`_.
+
 Use `examples/path_requests_run.py <examples/path_requests_run.py>`_ to run multiple optimizations as follows:
 
 .. code-block:: shell

docs/biblio.bib

@@ -874,7 +874,7 @@ month={Sept},}
   number = {7},
   journal = {Optics Express},
   urlyear = {2017-11-14},
-  year = {2012-03-26},
+  date = {2012-03-26},
   year = {2012},
   pages = {7777},
   author = {Bononi, A. and Serena, P. and Rossi, N. and Grellier, E. and Vacondio, F.}
@@ -1114,7 +1114,7 @@ month={Sept},}
   number = {26},
   journal = {Optics Express},
   urlyear = {2017-11-16},
-  year = {2013-12-30},
+  date = {2013-12-30},
   year = {2013},
   pages = {32254},
   author = {Bononi, Alberto and Beucher, Ottmar and Serena, Paolo}

docs/conf.py

@@ -173,5 +173,4 @@ texinfo_documents = [
      'Miscellaneous'),
 ]
 
-
+autodoc_default_flags = ['members', 'undoc-members', 'private-members', 'show-inheritance']
-

docs/source/gnpy.core.rst

@@ -4,10 +4,39 @@ gnpy\.core package
 Submodules
 ----------
 
+gnpy\.core\.ansi_escapes module
+-------------------------------
+
+.. automodule:: gnpy.core.ansi_escapes
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+gnpy\.core\.convert module
+--------------------------
+
+.. automodule:: gnpy.core.convert
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 gnpy\.core\.elements module
 ---------------------------
 
 .. automodule:: gnpy.core.elements
+
+gnpy\.core\.equipment module
+----------------------------
+
+.. automodule:: gnpy.core.equipment
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+gnpy\.core\.exceptions module
+-----------------------------
+
+.. automodule:: gnpy.core.exceptions
     :members:
     :undoc-members:
     :show-inheritance:
@@ -16,30 +45,34 @@ gnpy\.core\.execute module
 --------------------------
 
 .. automodule:: gnpy.core.execute
-    :members:
-    :undoc-members:
-    :show-inheritance:
 
 gnpy\.core\.info module
 -----------------------
 
 .. automodule:: gnpy.core.info
-    :members:
-    :undoc-members:
-    :show-inheritance:
 
 gnpy\.core\.network module
 --------------------------
 
 .. automodule:: gnpy.core.network
-    :members:
-    :undoc-members:
-    :show-inheritance:
 
 gnpy\.core\.node module
 -----------------------
 
 .. automodule:: gnpy.core.node
+
+gnpy\.core\.request module
+--------------------------
+
+.. automodule:: gnpy.core.request
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+gnpy\.core\.service_sheet module
+--------------------------------
+
+.. automodule:: gnpy.core.service_sheet
     :members:
     :undoc-members:
     :show-inheritance:
@@ -48,23 +81,14 @@ gnpy\.core\.units module
 ------------------------
 
 .. automodule:: gnpy.core.units
-    :members:
-    :undoc-members:
-    :show-inheritance:
 
 gnpy\.core\.utils module
 ------------------------
 
 .. automodule:: gnpy.core.utils
-    :members:
-    :undoc-members:
-    :show-inheritance:
 
 
 Module contents
 ---------------
 
 .. automodule:: gnpy.core
-    :members:
-    :undoc-members:
-    :show-inheritance:

docs/source/gnpy.rst

@@ -12,6 +12,3 @@ Module contents
 ---------------
 
 .. automodule:: gnpy
-    :members:
-    :undoc-members:
-    :show-inheritance:

examples/create_eqpt_sheet.py

@@ -11,64 +11,72 @@ If not present in the "Nodes" sheet, the "Type" column will be implicitly
 determined based on the topology.
 """
 
-from sys import exit
 try:
     from xlrd import open_workbook
 except ModuleNotFoundError:
     exit('Required: `pip install xlrd`')
 from argparse import ArgumentParser
-from collections import namedtuple, defaultdict
 
+PARSER = ArgumentParser()
+PARSER.add_argument('workbook', nargs='?', default='meshTopologyExampleV2.xls',
+                    help='create the mandatory columns in Eqpt sheet')
+ALL_ROWS = lambda sh, start=0: (sh.row(x) for x in range(start, sh.nrows))
 
-Shortlink = namedtuple('Link', 'src dest')
+class Node:
+    """ Node element contains uid, list of connected nodes and eqpt type
+    """
+    def __init__(self, uid, to_node):
+        self.uid = uid
+        self.to_node = to_node
+        self.eqpt = None
 
-Shortnode = namedtuple('Node', 'nodename eqt')
+    def __repr__(self):
+        return f'uid {self.uid} \nto_node {[node for node in self.to_node]}\neqpt {self.eqpt}\n'
 
-parser = ArgumentParser()
+    def __str__(self):
-parser.add_argument('workbook', nargs='?', default='meshTopologyExampleV2.xls',
+        return f'uid {self.uid} \nto_node {[node for node in self.to_node]}\neqpt {self.eqpt}\n'
-    help = 'create the mandatory columns in Eqpt sheet  ')
-all_rows = lambda sh, start=0: (sh.row(x) for x in range(start, sh.nrows))
 
 def read_excel(input_filename):
-    with open_workbook(input_filename) as wb:
+    """ read excel Nodes and Links sheets and create a dict of nodes with
+    their to_nodes and type of eqpt
+    """
+    with open_workbook(input_filename) as wobo:
         # reading Links sheet
-        links_sheet = wb.sheet_by_name('Links')
+        links_sheet = wobo.sheet_by_name('Links')
-        links = []
+        nodes = {}
-        nodeoccuranceinlinks = []
+        for row in ALL_ROWS(links_sheet, start=5):
-        links_by_src = defaultdict(list)
+            try:
-        links_by_dest = defaultdict(list)
+                nodes[row[0].value].to_node.append(row[1].value)
-        for row in all_rows(links_sheet, start=5):
+            except KeyError:
-            links.append(Shortlink(row[0].value,row[1].value))
+                nodes[row[0].value] = Node(row[0].value, [row[1].value])
-            links_by_src[row[0].value].append(Shortnode(row[1].value,''))
+            try:
-            links_by_dest[row[1].value].append(Shortnode(row[0].value,''))
+                nodes[row[1].value].to_node.append(row[0].value)
-            #print(f'source {links[len(links)-1].src} dest {links[len(links)-1].dest}')
+            except KeyError:
-            nodeoccuranceinlinks.append(row[0].value)
+                nodes[row[1].value] = Node(row[1].value, [row[0].value])
-            nodeoccuranceinlinks.append(row[1].value)
 
-        # reading Nodes sheet
+        nodes_sheet = wobo.sheet_by_name('Nodes')
-        nodes_sheet = wb.sheet_by_name('Nodes')
+        for row in ALL_ROWS(nodes_sheet, start=5):
-        nodes = []
+            node = row[0].value
-        node_degree = []
+            eqpt = row[6].value
-        for row in all_rows(nodes_sheet, start=5) :
+            try:
+                if eqpt == 'ILA' and len(nodes[node].to_node) != 2:
+                    print(f'Inconsistancy ILA node with degree > 2: {node} ')
+                    exit()
+                if eqpt == '' and len(nodes[node].to_node) == 2:
+                    nodes[node].eqpt = 'ILA'
+                elif eqpt == '' and len(nodes[node].to_node) != 2:
+                    nodes[node].eqpt = 'ROADM'
+                else:
+                    nodes[node].eqpt = eqpt
+            except KeyError:
+                print(f'inconsistancy between nodes and links sheet: {node} is not listed in links')
+                exit()
+        return nodes
 
-            temp_eqt = row[6].value
+def create_eqt_template(nodes, input_filename):
-            # verify node degree to confirm eqt type
+    """ writes list of node A node Z corresponding to Nodes and Links sheets in order
-            node_degree.append(nodeoccuranceinlinks.count(row[0].value))
+    to help user populating Eqpt
-            if temp_eqt.lower() == 'ila' and nodeoccuranceinlinks.count(row[0].value) !=2 :
+    """
-                print(f'Inconsistancy: node {nodes[len(nodes)-1]} has degree \
-                    {node_degree[len(nodes)-1]} and can not be an ILA ... replaced by ROADM')
-                temp_eqt = 'ROADM'
-            if temp_eqt == '' and nodeoccuranceinlinks.count(row[0].value) == 2 :
-                temp_eqt = 'ILA'
-            if temp_eqt == '' and nodeoccuranceinlinks.count(row[0].value) != 2 :
-                temp_eqt = 'ROADM'
-            # print(f'node {nodes[len(nodes)-1]} eqt {temp_eqt}')
-            nodes.append(Shortnode(row[0].value,temp_eqt))
-            # print(len(nodes)-1)
-            print(f'reading: node {nodes[len(nodes)-1].nodename} eqpt {temp_eqt}')
-        return links,nodes, links_by_src , links_by_dest
-
-def create_eqt_template(links,nodes, links_by_src , links_by_dest, input_filename):
     output_filename = f'{input_filename[:-4]}_eqpt_sheet.txt'
     with open(output_filename, 'w', encoding='utf-8') as my_file:
         # print header similar to excel
@@ -77,27 +85,17 @@ def create_eqt_template(links,nodes, links_by_src , links_by_dest, input_filenam
            \nNode A \tNode Z \tamp type \tatt_in \tamp gain \ttilt \tatt_out\
            amp type   \tatt_in \tamp gain   \ttilt   \tatt_out\n')
 
-        tab = []
+
-        temp = []
+        for node in nodes.values():
-        i = 0
+            if node.eqpt == 'ILA':
-        for lk in links:
+                my_file.write(f'{node.uid}\t{node.to_node[0]}\n')
-            if [e for n,e in nodes if n==lk.src][0] != 'FUSED' :
+            if node.eqpt == 'ROADM':
-                temp = [lk.src , lk.dest]
+                for to_node in node.to_node:
-                tab.append(temp)
+                    my_file.write(f'{node.uid}\t{to_node}\n')
-                my_file.write(f'{temp[0]}\t{temp[1]}\n')
+
-        for n in nodes :
-            if n.eqt.lower() == 'roadm' :
-                for src in  links_by_dest[n.nodename] :
-                    temp = [n.nodename , src.nodename]
-                    tab.append(temp)
-                    # print(temp)
-                    my_file.write(f'{temp[0]}\t{temp[1]}\n')
-            i = i + 1
         print(f'File {output_filename} successfully created with Node A - Node Z ' +
-        ' entries for Eqpt sheet in excel file.')
+              ' entries for Eqpt sheet in excel file.')
 
 if __name__ == '__main__':
-    args = parser.parse_args()
+    ARGS = PARSER.parse_args()
-    input_filename = args.workbook
+    create_eqt_template(read_excel(ARGS.workbook), ARGS.workbook)
-    links,nodes,links_by_src, links_by_dest = read_excel(input_filename)
-    create_eqt_template(links,nodes, links_by_src , links_by_dest , input_filename)

examples/eqpt_config.json

@@ -159,6 +159,36 @@
             "gamma": 0.000843
             }
       ],
+      "RamanFiber":[{
+            "type_variety": "SSMF",
+            "dispersion": 1.67e-05,
+            "gamma": 0.00127,
+            "raman_efficiency": {
+              "cr":[
+                  0, 9.4E-06, 2.92E-05, 4.88E-05, 6.82E-05, 8.31E-05, 9.4E-05, 0.0001014, 0.0001069, 0.0001119,
+                  0.0001217, 0.0001268, 0.0001365, 0.000149, 0.000165, 0.000181, 0.0001977, 0.0002192, 0.0002469,
+                  0.0002749, 0.0002999, 0.0003206, 0.0003405, 0.0003592, 0.000374, 0.0003826, 0.0003841, 0.0003826,
+                  0.0003802, 0.0003756, 0.0003549, 0.0003795, 0.000344, 0.0002933, 0.0002024, 0.0001158, 8.46E-05,
+                  7.14E-05, 6.86E-05, 8.5E-05, 8.93E-05, 9.01E-05, 8.15E-05, 6.67E-05, 4.37E-05, 3.28E-05, 2.96E-05,
+                  2.65E-05, 2.57E-05, 2.81E-05, 3.08E-05, 3.67E-05, 5.85E-05, 6.63E-05, 6.36E-05, 5.5E-05, 4.06E-05,
+                  2.77E-05, 2.42E-05, 1.87E-05, 1.6E-05, 1.4E-05, 1.13E-05, 1.05E-05, 9.8E-06, 9.8E-06, 1.13E-05,
+                  1.64E-05, 1.95E-05, 2.38E-05, 2.26E-05, 2.03E-05, 1.48E-05, 1.09E-05, 9.8E-06, 1.05E-05, 1.17E-05,
+                  1.25E-05, 1.21E-05, 1.09E-05, 9.8E-06, 8.2E-06, 6.6E-06, 4.7E-06, 2.7E-06, 1.9E-06, 1.2E-06, 4E-07,
+                  2E-07, 1E-07
+              ],
+              "frequency_offset":[
+                0, 0.5e12, 1e12, 1.5e12, 2e12, 2.5e12, 3e12, 3.5e12, 4e12, 4.5e12, 5e12, 5.5e12, 6e12, 6.5e12, 7e12,
+                7.5e12, 8e12, 8.5e12, 9e12, 9.5e12, 10e12, 10.5e12, 11e12, 11.5e12, 12e12, 12.5e12, 12.75e12,
+                13e12, 13.25e12, 13.5e12, 14e12, 14.5e12, 14.75e12, 15e12, 15.5e12, 16e12, 16.5e12, 17e12,
+                17.5e12, 18e12, 18.25e12, 18.5e12, 18.75e12, 19e12, 19.5e12, 20e12, 20.5e12, 21e12, 21.5e12,
+                22e12, 22.5e12, 23e12, 23.5e12, 24e12, 24.5e12, 25e12, 25.5e12, 26e12, 26.5e12, 27e12, 27.5e12, 28e12,
+                28.5e12, 29e12, 29.5e12, 30e12, 30.5e12, 31e12, 31.5e12, 32e12, 32.5e12, 33e12, 33.5e12, 34e12, 34.5e12,
+                35e12, 35.5e12, 36e12, 36.5e12, 37e12, 37.5e12, 38e12, 38.5e12, 39e12, 39.5e12, 40e12, 40.5e12, 41e12,
+                41.5e12, 42e12
+              ]
+              }
+            }
+      ],
       "Span":[{
             "power_mode":true,
             "delta_power_range_db": [-2,3,0.5],
@@ -177,8 +207,8 @@
             "target_pch_out_db": -20,
             "add_drop_osnr": 38,
             "restrictions": {
-                            "preamp_variety_list":["low_gain_preamp", "high_gain_preamp"],
+                            "preamp_variety_list":[],
-                            "booster_variety_list":["std_booster"]
+                            "booster_variety_list":[]
                             }            
             }],
       "SI":[{

examples/meshTopologyExampleV2.json

@@ -681,25 +681,6 @@
         "out_voa": null
       }
     },
-    {
-      "uid": "east edfa in Lorient_KMA to Vannes_KBE",
-      "metadata": {
-        "location": {
-          "city": "Lorient_KMA",
-          "region": "RLD",
-          "latitude": 2.0,
-          "longitude": 3.0
-        }
-      },
-      "type": "Edfa",
-      "type_variety": "std_low_gain",
-      "operational": {
-        "gain_target": null,
-        "delta_p": 1.0,
-        "tilt_target": 0,
-        "out_voa": null
-      }
-    },
     {
       "uid": "east edfa in Lannion_CAS to Stbrieuc",
       "metadata": {
@@ -1041,6 +1022,21 @@
         "tilt_target": 0,
         "out_voa": null
       }
+    },
+    {
+      "uid": "east edfa in Lorient_KMA to Vannes_KBE",
+      "metadata": {
+        "location": {
+          "city": "Lorient_KMA",
+          "region": "RLD",
+          "latitude": 2.0,
+          "longitude": 3.0
+        }
+      },
+      "type": "Fused",
+      "params": {
+        "loss": 0
+      }
     }
   ],
   "connections": [

examples/path_requests_run.py

@@ -30,10 +30,11 @@ from gnpy.core.utils import db2lin, lin2db
 from gnpy.core.request import (Path_request, Result_element, compute_constrained_path,
                               propagate, jsontocsv, Disjunction, compute_path_dsjctn, requests_aggregation,
                               propagate_and_optimize_mode)
+from gnpy.core.exceptions import ConfigurationError, EquipmentConfigError, NetworkTopologyError
+import gnpy.core.ansi_escapes as ansi_escapes
 from copy import copy, deepcopy
 from textwrap import dedent
 from math import ceil
-import time
 
 #EQPT_LIBRARY_FILENAME = Path(__file__).parent / 'eqpt_config.json'
 
@@ -142,44 +143,6 @@ def load_requests(filename,eqpt_filename):
             json_data = loads(f.read())
     return json_data
 
-def compute_path(network, equipment, pathreqlist):
-
-    # This function is obsolete and not relevant with respect to network building: suggest either to correct
-    # or to suppress it
-    
-    path_res_list = []
-
-    for pathreq in pathreqlist:
-        #need to rebuid the network for each path because the total power
-        #can be different and the choice of amplifiers in autodesign is power dependant
-        #but the design is the same if the total power is the same
-        #TODO parametrize the total spectrum power so the same design can be shared
-        p_db = lin2db(pathreq.power*1e3)
-        p_total_db = p_db + lin2db(pathreq.nb_channel)
-        build_network(network, equipment, p_db, p_total_db)
-        pathreq.nodes_list.append(pathreq.destination)
-        #we assume that the destination is a strict constraint
-        pathreq.loose_list.append('strict')
-        print(f'Computing path from {pathreq.source} to {pathreq.destination}')
-        print(f'with path constraint: {[pathreq.source]+pathreq.nodes_list}') #adding first node to be clearer on the output
-        total_path = compute_constrained_path(network, pathreq)
-        print(f'Computed path (roadms):{[e.uid for e in total_path  if isinstance(e, Roadm)]}\n')
-
-        if total_path :
-            total_path = propagate(total_path,pathreq,equipment, show=False)
-        else:
-            total_path = []
-        # we record the last tranceiver object in order to have th whole
-        # information about spectrum. Important Note: since transceivers
-        # attached to roadms are actually logical elements to simulate
-        # performance, several demands having the same destination may use
-        # the same transponder for the performance simaulation. This is why
-        # we use deepcopy: to ensure each propagation is recorded and not
-        # overwritten
-
-        path_res_list.append(deepcopy(total_path))
-    return path_res_list
-
 def compute_path_with_disjunction(network, equipment, pathreqlist, pathlist):
     
     # use a list but a dictionnary might be helpful to find path bathsed on request_id
@@ -203,7 +166,8 @@ def compute_path_with_disjunction(network, equipment, pathreqlist, pathlist):
         # print(f'{pathreq.baud_rate}   {pathreq.power}   {pathreq.spacing}   {pathreq.nb_channel}')
         if total_path :
             if pathreq.baud_rate is not None:
-                total_path = propagate(total_path,pathreq,equipment, show=False)
+                total_path = propagate(total_path,pathreq,equipment)
+                # for el in total_path: print(el)
                 temp_snr01nm = round(mean(total_path[-1].snr+lin2db(pathreq.baud_rate/(12.5e9))),2)
                 if temp_snr01nm < pathreq.OSNR :
                     msg = f'\tWarning! Request {pathreq.request_id} computed path from {pathreq.source} to {pathreq.destination} does not pass with {pathreq.tsp_mode}\n' +\
@@ -299,22 +263,31 @@ def path_result_json(pathresult):
 
 
 if __name__ == '__main__':
-    start = time.time()
     args = parser.parse_args()
     basicConfig(level={2: DEBUG, 1: INFO, 0: CRITICAL}.get(args.verbose, DEBUG))
     logger.info(f'Computing path requests {args.service_filename} into JSON format')
     print('\x1b[1;34;40m'+f'Computing path requests {args.service_filename} into JSON format'+ '\x1b[0m')
     # for debug
     # print( args.eqpt_filename)
-    data = load_requests(args.service_filename,args.eqpt_filename)
+    try:
-    equipment = load_equipment(args.eqpt_filename)
+        data = load_requests(args.service_filename,args.eqpt_filename)
-    network = load_network(args.network_filename,equipment)
+        equipment = load_equipment(args.eqpt_filename)
+        network = load_network(args.network_filename,equipment)
+    except EquipmentConfigError as e:
+        print(f'{ansi_escapes.red}Configuration error in the equipment library:{ansi_escapes.reset} {e}')
+        exit(1)
+    except NetworkTopologyError as e:
+        print(f'{ansi_escapes.red}Invalid network definition:{ansi_escapes.reset} {e}')
+        exit(1)
+    except ConfigurationError as e:
+        print(f'{ansi_escapes.red}Configuration error:{ansi_escapes.reset} {e}')
+        exit(1)
 
     # Build the network once using the default power defined in SI in eqpt config
     # TODO power density : db2linp(ower_dbm": 0)/power_dbm": 0 * nb channels as defined by
-    # spacing, f_min and f_max 
+    # spacing, f_min and f_max
     p_db = equipment['SI']['default'].power_dbm
-    
+
     p_total_db = p_db + lin2db(automatic_nch(equipment['SI']['default'].f_min,\
         equipment['SI']['default'].f_max, equipment['SI']['default'].spacing))
     build_network(network, equipment, p_db, p_total_db)
@@ -322,7 +295,7 @@ if __name__ == '__main__':
 
     rqs = requests_from_json(data, equipment)
 
-    # check that request ids are unique. Non unique ids, may 
+    # check that request ids are unique. Non unique ids, may
     # mess the computation : better to stop the computation
     all_ids = [r.request_id for r in rqs]
     if len(all_ids) != len(set(all_ids)):
@@ -341,7 +314,7 @@ if __name__ == '__main__':
     # need to warn or correct in case of wrong disjunction form
     # disjunction must not be repeated with same or different ids
     dsjn = correct_disjn(dsjn)
-        
+
     # Aggregate demands with same exact constraints
     print('\x1b[1;34;40m'+f'Aggregating similar requests'+ '\x1b[0m')
 
@@ -350,17 +323,15 @@ if __name__ == '__main__':
 
     print('\x1b[1;34;40m'+'The following services have been requested:'+ '\x1b[0m')
     print(rqs)
-    
+
     print('\x1b[1;34;40m'+f'Computing all paths with constraints'+ '\x1b[0m')
     pths = compute_path_dsjctn(network, equipment, rqs, dsjn)
 
     print('\x1b[1;34;40m'+f'Propagating on selected path'+ '\x1b[0m')
     propagatedpths = compute_path_with_disjunction(network, equipment, rqs, pths)
 
-    end = time.time()
-    print(f'computation time {end-start}')
     print('\x1b[1;34;40m'+f'Result summary'+ '\x1b[0m')
-    
+
     header = ['req id', '  demand','  snr@bandwidth','  snr@0.1nm','  Receiver minOSNR', '  mode', '  Gbit/s' , '  nb of tsp pairs']
     data = []
     data.append(header)
@@ -377,7 +348,7 @@ if __name__ == '__main__':
     firstcol_width = max(len(row[0]) for row in data )   # padding
     secondcol_width = max(len(row[1]) for row in data )   # padding
     for row in data:
-        firstcol = ''.join(row[0].ljust(firstcol_width)) 
+        firstcol = ''.join(row[0].ljust(firstcol_width))
         secondcol = ''.join(row[1].ljust(secondcol_width))
         remainingcols = ''.join(word.center(col_width,' ') for word in row[2:])
         print(f'{firstcol} {secondcol} {remainingcols}')
@@ -396,4 +367,3 @@ if __name__ == '__main__':
             with open(fnamecsv,"w", encoding='utf-8') as fcsv :
                 jsontocsv(temp,equipment,fcsv)
                 print('\x1b[1;34;40m'+f'saving in {args.output} and {fnamecsv}'+ '\x1b[0m')
-

examples/raman_edfa_example_network.json

@@ -0,0 +1,98 @@
+{
+  "elements": [
+    {
+      "uid": "Site_A",
+      "type": "Transceiver",
+      "metadata": {
+        "location": {
+          "latitude": 0,
+          "longitude": 0,
+          "city": "Site A",
+          "region": ""
+        }
+      }
+    },
+    {
+      "uid": "Span1",
+      "type": "RamanFiber",
+      "type_variety": "SSMF",
+      "operational": {
+        "temperature": 283,
+        "raman_pumps": [
+          {
+            "power": 200e-3,
+            "frequency": 205e12,
+            "propagation_direction": "counterprop"
+          },
+          {
+            "power": 206e-3,
+            "frequency": 201e12,
+            "propagation_direction": "counterprop"
+          }
+        ]
+      },
+      "params": {
+        "type_variety": "SSMF",
+        "length": 80.0,
+        "loss_coef": 0.2,
+        "length_units": "km",
+        "att_in": 0,
+        "con_in": 0.5,
+        "con_out": 0.5
+      },
+      "metadata": {
+        "location": {
+          "latitude": 1,
+          "longitude": 0,
+          "city": null,
+          "region": ""
+        }
+      }
+    },
+    {
+      "uid": "Edfa1",
+      "type": "Edfa",
+      "type_variety": "std_low_gain",
+      "operational": {
+        "gain_target": 15.0,
+        "delta_p": -2,
+        "tilt_target": 0,
+        "out_voa": 0
+      },
+      "metadata": {
+        "location": {
+          "latitude": 2,
+          "longitude": 0,
+          "city": null,
+          "region": ""
+        }
+      }
+    },
+    {
+      "uid": "Site_B",
+      "type": "Transceiver",
+      "metadata": {
+        "location": {
+          "latitude": 2,
+          "longitude": 0,
+          "city": "Site B",
+          "region": ""
+        }
+      }
+    }
+  ],
+  "connections": [
+    {
+      "from_node": "Site_A",
+      "to_node": "Span1"
+    },
+    {
+      "from_node": "Span1",
+      "to_node": "Edfa1"
+    },
+    {
+      "from_node": "Edfa1",
+      "to_node": "Site_B"
+    }
+  ]
+}

examples/sim_params.json

@@ -0,0 +1,14 @@
+{
+  "raman_computed_channels": [1, 18, 37, 56, 75],
+  "raman_parameters": {
+    "flag_raman": true,
+    "space_resolution": 10e3,
+    "tolerance": 1e-8
+  },
+  "nli_parameters": {
+  	"nli_method_name": "ggn_spectrally_separated",
+  	"wdm_grid_size": 50e9,
+  	"dispersion_tolerance": 1,
+  	"phase_shift_tollerance": 0.1
+  }
+}

examples/transmission_main_example.py

@@ -18,14 +18,16 @@ from pathlib import Path
 from json import loads
 from collections import Counter
 from logging import getLogger, basicConfig, INFO, ERROR, DEBUG
-from numpy import linspace, mean
+from numpy import linspace, mean, log10
 from matplotlib.pyplot import show, axis, figure, title, text
 from networkx import (draw_networkx_nodes, draw_networkx_edges,
                       draw_networkx_labels, dijkstra_path)
-from gnpy.core.network import load_network, build_network, save_network
+from gnpy.core.network import load_network, build_network, save_network, load_sim_params, configure_network
-from gnpy.core.elements import Transceiver, Fiber, Edfa, Roadm
+from gnpy.core.elements import Transceiver, Fiber, RamanFiber, Edfa, Roadm
 from gnpy.core.info import create_input_spectral_information, SpectralInformation, Channel, Power, Pref
 from gnpy.core.request import Path_request, RequestParams, compute_constrained_path, propagate2
+from gnpy.core.exceptions import ConfigurationError, EquipmentConfigError, NetworkTopologyError
+import gnpy.core.ansi_escapes as ansi_escapes
 
 logger = getLogger(__name__)
 
@@ -97,7 +99,7 @@ def plot_results(network, path, source, destination, infos):
     show()
 
 
-def main(network, equipment, source, destination, req = None):
+def main(network, equipment, source, destination, sim_params, req=None):
     result_dicts = {}
     network_data = [{
                     'network_name'  : str(args.filename),
@@ -124,8 +126,14 @@ def main(network, equipment, source, destination, req = None):
     build_network(network, equipment, pref_ch_db, pref_total_db)
     path = compute_constrained_path(network, req)
 
-    spans = [s.length for s in path if isinstance(s, Fiber)]
+    if len([s.length for s in path if isinstance(s, RamanFiber)]):
-    print(f'\nThere are {len(spans)} fiber spans over {sum(spans):.0f}m between {source.uid} and {destination.uid}')
+        if sim_params is None:
+            print(f'{ansi_escapes.red}Invocation error:{ansi_escapes.reset} RamanFiber requires passing simulation params via --sim-params')
+            exit(1)
+        configure_network(network, sim_params)
+
+    spans = [s.length for s in path if isinstance(s, RamanFiber) or isinstance(s, Fiber)]
+    print(f'\nThere are {len(spans)} fiber spans over {sum(spans)/1000:.0f} km between {source.uid} and {destination.uid}')
     print(f'\nNow propagating between {source.uid} and {destination.uid}:')
 
     try:
@@ -142,16 +150,18 @@ def main(network, equipment, source, destination, req = None):
     for dp_db in power_range:
         req.power = db2lin(pref_ch_db + dp_db)*1e-3
         if power_mode:
-            print(f'\nPropagating with input power = {lin2db(req.power*1e3):.2f}dBm :')
+            print(f'\nPropagating with input power = {ansi_escapes.cyan}{lin2db(req.power*1e3):.2f} dBm{ansi_escapes.reset}:')
         else:
-            print(f'\nPropagating in gain mode: power cannot be set manually')
+            print(f'\nPropagating in {ansi_escapes.cyan}gain mode{ansi_escapes.reset}: power cannot be set manually')
-        infos = propagate2(path, req, equipment, show=len(power_range)==1)
+        infos = propagate2(path, req, equipment)
+        if len(power_range) == 1:
+            for elem in path:
+                print(elem)
         if power_mode:
-            print(f'\nTransmission result for input power = {lin2db(req.power*1e3):.2f}dBm :')
+            print(f'\nTransmission result for input power = {lin2db(req.power*1e3):.2f} dBm:')
         else:
             print(f'\nTransmission results:')
-        #info message in gain mode
+        print(f'  Final SNR total (signal bw): {ansi_escapes.cyan}{mean(destination.snr):.02f} dB{ansi_escapes.reset}')
-        print(destination)
 
         #print(f'\n !!!!!!!!!!!!!!!!!     TEST POINT         !!!!!!!!!!!!!!!!!!!!!')
         #print(f'carriers ase output of {path[1]} =\n {list(path[1].carriers("out", "nli"))}')
@@ -172,8 +182,7 @@ def main(network, equipment, source, destination, req = None):
                         'OSNR_ASE_signal_bw'    : round(mean(destination.osnr_ase),2),
                         'SNR_nli_signal_bw'     : round(mean(destination.osnr_nli),2),
                         'SNR_total_signal_bw'   : round(mean(destination.snr),2)
-                                })          
+                                })
-        #info message in gain mode
     write_csv(result_dicts, 'simulation_result.csv')
     return path, infos
 
@@ -181,6 +190,9 @@ def main(network, equipment, source, destination, req = None):
 parser = ArgumentParser()
 parser.add_argument('-e', '--equipment', type=Path,
                     default=Path(__file__).parent / 'eqpt_config.json')
+parser.add_argument('--sim-params', type=Path,
+                    default=None, help='Path to the JSON containing simulation parameters (required for Raman)')
+parser.add_argument('--show-channels', action='store_true', help='Show final per-channel OSNR summary')
 parser.add_argument('-pl', '--plot', action='store_true')
 parser.add_argument('-v', '--verbose', action='count', default=0, help='increases verbosity for each occurence')
 parser.add_argument('-l', '--list-nodes', action='store_true', help='list all transceiver nodes')
@@ -196,8 +208,19 @@ if __name__ == '__main__':
     args = parser.parse_args()
     basicConfig(level={0: ERROR, 1: INFO, 2: DEBUG}.get(args.verbose, DEBUG))
 
-    equipment = load_equipment(args.equipment)
+    try:
-    network = load_network(args.filename, equipment, args.names_matching)
+        equipment = load_equipment(args.equipment)
+        network = load_network(args.filename, equipment, args.names_matching)
+        sim_params = load_sim_params(args.sim_params) if args.sim_params is not None else None
+    except EquipmentConfigError as e:
+        print(f'{ansi_escapes.red}Configuration error in the equipment library:{ansi_escapes.reset} {e}')
+        exit(1)
+    except NetworkTopologyError as e:
+        print(f'{ansi_escapes.red}Invalid network definition:{ansi_escapes.reset} {e}')
+        exit(1)
+    except ConfigurationError as e:
+        print(f'{ansi_escapes.red}Configuration error:{ansi_escapes.reset} {e}')
+        exit(1)
 
     if args.plot:
         plot_baseline(network)
@@ -266,9 +289,19 @@ if __name__ == '__main__':
         trx_params['power'] = db2lin(float(args.power))*1e-3
     params.update(trx_params)
     req = Path_request(**params)
-    path, infos = main(network, equipment, source, destination, req)
+    path, infos = main(network, equipment, source, destination, sim_params, req)
     save_network(args.filename, network)
 
+    if args.show_channels:
+        print('\nThe total SNR per channel at the end of the line is:')
+        print('{:>5}{:>26}{:>26}{:>28}{:>28}{:>28}' \
+            .format('Ch. #', 'Channel frequency (THz)', 'Channel power (dBm)', 'OSNR ASE (signal bw, dB)', 'SNR NLI (signal bw, dB)', 'SNR total (signal bw, dB)'))
+        for final_carrier, ch_osnr, ch_snr_nl, ch_snr in zip(infos[path[-1]][1].carriers, path[-1].osnr_ase, path[-1].osnr_nli, path[-1].snr):
+            ch_freq = final_carrier.frequency * 1e-12
+            ch_power = lin2db(final_carrier.power.signal*1e3)
+            print('{:5}{:26.2f}{:26.2f}{:28.2f}{:28.2f}{:28.2f}' \
+                .format(final_carrier.channel_number, round(ch_freq, 2), round(ch_power, 2), round(ch_osnr, 2), round(ch_snr_nl, 2), round(ch_snr, 2)))
+
     if not args.source:
         print(f'\n(No source node specified: picked {source.uid})')
     elif not valid_source:

gnpy/core/ansi_escapes.py

@@ -0,0 +1,13 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+'''
+gnpy.core.ansi_escapes
+======================
+
+A random subset of ANSI terminal escape codes for colored messages
+'''
+
+red = '\x1b[1;31;40m'
+cyan = '\x1b[1;36;40m'
+reset = '\x1b[0m'

gnpy/core/convert.py

@@ -31,6 +31,7 @@ from itertools import chain
 from json import dumps
 from pathlib import Path
 from difflib import get_close_matches
+from gnpy.core.utils import silent_remove
 import time
 
 all_rows = lambda sh, start=0: (sh.row(x) for x in range(start, sh.nrows))
@@ -54,7 +55,9 @@ class Node(object):
         'region':       '',
         'latitude':     0,
         'longitude':    0,
-        'node_type':    'ILA'
+        'node_type':    'ILA',
+        'booster_restriction' : '',
+        'preamp_restriction'  : ''
     }
 
 class Link(object):
@@ -235,7 +238,6 @@ def sanity_check(nodes, links, nodes_by_city, links_by_city, eqpts_by_city):
 
 def convert_file(input_filename, names_matching=False, filter_region=[]):
     nodes, links, eqpts = parse_excel(input_filename)
-
     if filter_region:
         nodes = [n for n in nodes if n.region.lower() in filter_region]
         cities = {n.city for n in nodes}
@@ -244,10 +246,8 @@ def convert_file(input_filename, names_matching=False, filter_region=[]):
         cities = {lnk.from_city for lnk in links} | {lnk.to_city for lnk in links}
         nodes = [n for n in nodes if n.city in cities]
 
-
     global nodes_by_city
     nodes_by_city = {n.city: n for n in nodes}
-
     #create matching dictionary for node name mismatch analysis
 
     cities = {''.join(c.strip() for c in n.city.split('C+L')).lower(): n.city for n in nodes}
@@ -298,7 +298,22 @@ def convert_file(input_filename, names_matching=False, filter_region=[]):
                                         'latitude':  x.latitude,
                                         'longitude': x.longitude}},
               'type': 'Roadm'}
-             for x in nodes_by_city.values() if x.node_type.lower() == 'roadm'] +
+             for x in nodes_by_city.values() if x.node_type.lower() == 'roadm' \
+                 and x.booster_restriction == '' and x.preamp_restriction == ''] +
+            [{'uid': f'roadm {x.city}',
+              'params' : {
+                'restrictions': {
+                  'preamp_variety_list': silent_remove(x.preamp_restriction.split(' | '),''),
+                  'booster_variety_list': silent_remove(x.booster_restriction.split(' | '),'')
+                  }
+              },
+              'metadata': {'location': {'city':      x.city,
+                                        'region':    x.region,
+                                        'latitude':  x.latitude,
+                                        'longitude': x.longitude}},
+              'type': 'Roadm'}
+             for x in nodes_by_city.values() if x.node_type.lower() == 'roadm' and \
+                 (x.booster_restriction != '' or x.preamp_restriction != '')] +
             [{'uid': f'west fused spans in {x.city}',
               'metadata': {'location': {'city':      x.city,
                                         'region':    x.region,
@@ -348,8 +363,9 @@ def convert_file(input_filename, names_matching=False, filter_region=[]):
                               'delta_p':     e.east_amp_dp,
                               'tilt_target': e.east_tilt,
                               'out_voa'    : e.east_att_out}
-            }
+             }
-             for e in eqpts if e.east_amp_type.lower() != ''] +
+             for e in eqpts if (e.east_amp_type.lower() != '' and \
+                                e.east_amp_type.lower() != 'fused')] +
             [{'uid': f'west edfa in {e.from_city} to {e.to_city}',
               'metadata': {'location': {'city':      nodes_by_city[e.from_city].city,
                                         'region':    nodes_by_city[e.from_city].region,
@@ -361,8 +377,31 @@ def convert_file(input_filename, names_matching=False, filter_region=[]):
                               'delta_p':     e.west_amp_dp,
                               'tilt_target': e.west_tilt,
                               'out_voa'    : e.west_att_out}
-              }
+             }
-             for e in eqpts if e.west_amp_type.lower() != ''],
+             for e in eqpts if (e.west_amp_type.lower() != '' and \
+                                e.west_amp_type.lower() != 'fused')] +
+            # fused edfa variety is a hack to indicate that there should not be
+            # booster amplifier out the roadm.
+            # If user specifies ILA in Nodes sheet and fused in Eqpt sheet, then assumes that
+            # this is a fused nodes.
+            [{'uid': f'east edfa in {e.from_city} to {e.to_city}',
+              'metadata': {'location': {'city':      nodes_by_city[e.from_city].city,
+                                        'region':    nodes_by_city[e.from_city].region,
+                                        'latitude':  nodes_by_city[e.from_city].latitude,
+                                        'longitude': nodes_by_city[e.from_city].longitude}},
+              'type': 'Fused',
+              'params': {'loss': 0}
+             }
+             for e in eqpts if e.east_amp_type.lower() == 'fused'] +
+            [{'uid': f'west edfa in {e.from_city} to {e.to_city}',
+              'metadata': {'location': {'city':      nodes_by_city[e.from_city].city,
+                                        'region':    nodes_by_city[e.from_city].region,
+                                        'latitude':  nodes_by_city[e.from_city].latitude,
+                                        'longitude': nodes_by_city[e.from_city].longitude}},
+              'type': 'Fused',
+              'params': {'loss': 0}
+             }
+             for e in eqpts if e.west_amp_type.lower() == 'fused'],
         'connections':
             list(chain.from_iterable([eqpt_connection_by_city(n.city)
             for n in nodes]))
@@ -414,7 +453,9 @@ def parse_excel(input_filename):
         'Region':       'region',
         'Latitude':     'latitude',
         'Longitude':    'longitude',
-        'Type':         'node_type'
+        'Type':         'node_type',
+        'Booster_restriction': 'booster_restriction',
+        'Preamp_restriction': 'preamp_restriction'
     }
     eqpt_headers = \
     {  'Node A': 'from_city',
@@ -467,10 +508,10 @@ def parse_excel(input_filename):
     # sanity check
     all_cities = Counter(n.city for n in nodes)
     if len(all_cities) != len(nodes):
-        ValueError(f'Duplicate city: {all_cities}')
+        raise ValueError(f'Duplicate city: {all_cities}')
     if any(ln.from_city not in all_cities or
            ln.to_city   not in all_cities for ln in links):
-        ValueError(f'Bad link.')
+        raise ValueError(f'Bad link.')
 
     return nodes, links, eqpts
 
@@ -571,7 +612,7 @@ def midpoint(city_a, city_b):
 #output_json_file_name = 'coronet_conus_example.json'
 #TODO get column size automatically from tupple size
 
-NODES_COLUMN = 8
+NODES_COLUMN = 10
 NODES_LINE = 4
 LINKS_COLUMN = 16
 LINKS_LINE = 3

gnpy/core/elements.py

@@ -7,18 +7,18 @@ gnpy.core.elements
 
 This module contains standard network elements.
 
-A network element is a Python callable. It takes a .info.SpectralInformation
+A network element is a Python callable. It takes a :class:`.info.SpectralInformation`
 object and returns a copy with appropriate fields affected. This structure
 represents spectral information that is "propogated" by this network element.
 Network elements must have only a local "view" of the network and propogate
-SpectralInformation using only this information. They should be independent and
+:class:`.info.SpectralInformation` using only this information. They should be independent and
 self-contained.
 
 Network elements MUST implement two attributes .uid and .name representing a
 unique identifier and a printable name.
 '''
 
-from numpy import abs, arange, arcsinh, array, exp, divide, errstate
+from numpy import abs, arange, array, exp, divide, errstate
 from numpy import interp, log10, mean, pi, polyfit, polyval, sum
 from scipy.constants import c, h
 from collections import namedtuple
@@ -26,6 +26,7 @@ from collections import namedtuple
 from gnpy.core.node import Node
 from gnpy.core.units import UNITS
 from gnpy.core.utils import lin2db, db2lin, itufs, itufl, snr_sum
+from gnpy.core.science_utils import propagate_raman_fiber, _psi
 
 class Transceiver(Node):
     def __init__(self, *args, **kwargs):
@@ -117,7 +118,7 @@ class Transceiver(Node):
         self._calc_snr(spectral_info)
         return spectral_info
 
-RoadmParams = namedtuple('RoadmParams', 'target_pch_out_db add_drop_osnr')
+RoadmParams = namedtuple('RoadmParams', 'target_pch_out_db add_drop_osnr restrictions')
 
 class Roadm(Node):
     def __init__(self, *args, params, **kwargs):
@@ -126,15 +127,19 @@ class Roadm(Node):
         self.effective_loss = None
         self.effective_pch_out_db = self.params.target_pch_out_db
         self.passive = True
+        self.restrictions = self.params.restrictions
 
     @property
     def to_json(self):
         return {'uid'       : self.uid,
                 'type'      : type(self).__name__,
-                'params'    : {'target_pch_out_db' : self.effective_pch_out_db},
+                'params'    : {
+                    'target_pch_out_db' : self.effective_pch_out_db,
+                    'restrictions'      : self.restrictions
+                    },
                 'metadata'      : {
                     'location': self.metadata['location']._asdict()
-                                    }
+                                }
                 }
 
     def __repr__(self):
@@ -150,8 +155,8 @@ class Roadm(Node):
         #all ingress channels in xpress are set to this power level
         #but add channels are not, so we define an effective loss
         #in the case of add channels
-        self.effective_pch_out_db = min(pref.pi, self.params.target_pch_out_db)
+        self.effective_pch_out_db = min(pref.p_spani, self.params.target_pch_out_db)
-        self.effective_loss = pref.pi - self.effective_pch_out_db
+        self.effective_loss = pref.p_spani - self.effective_pch_out_db
         carriers_power = array([c.power.signal +c.power.nli+c.power.ase for c in carriers])
         carriers_att = list(map(lambda x : lin2db(x*1e3)-self.params.target_pch_out_db, carriers_power))
         exceeding_att = -min(list(filter(lambda x: x < 0, carriers_att)), default = 0)
@@ -159,17 +164,17 @@ class Roadm(Node):
         for carrier_att, carrier in zip(carriers_att, carriers) :
             pwr = carrier.power
             pwr = pwr._replace( signal = pwr.signal/carrier_att,
-                                nonlinear_interference = pwr.nli/carrier_att,
+                                nli = pwr.nli/carrier_att,
-                                amplified_spontaneous_emission = pwr.ase/carrier_att)
+                                ase = pwr.ase/carrier_att)
             yield carrier._replace(power=pwr)
 
     def update_pref(self, pref):
-        return pref._replace(p_span0=pref.p0, p_spani=self.effective_pch_out_db)
+        return pref._replace(p_span0=pref.p_span0, p_spani=self.effective_pch_out_db)
 
     def __call__(self, spectral_info):
         carriers = tuple(self.propagate(spectral_info.pref, *spectral_info.carriers))
         pref = self.update_pref(spectral_info.pref)
-        return spectral_info.update(carriers=carriers, pref=pref)
+        return spectral_info._replace(carriers=carriers, pref=pref)
 
 FusedParams = namedtuple('FusedParams', 'loss')
 
@@ -186,6 +191,9 @@ class Fused(Node):
     def to_json(self):
         return {'uid'       : self.uid,
                 'type'      : type(self).__name__,
+                'params'    :{
+                    'loss': self.loss
+                },
                 'metadata'      : {
                     'location': self.metadata['location']._asdict()
                                     }
@@ -204,17 +212,17 @@ class Fused(Node):
         for carrier in carriers:
             pwr = carrier.power
             pwr = pwr._replace(signal=pwr.signal/attenuation,
-                               nonlinear_interference=pwr.nli/attenuation,
+                               nli=pwr.nli/attenuation,
-                               amplified_spontaneous_emission=pwr.ase/attenuation)
+                               ase=pwr.ase/attenuation)
             yield carrier._replace(power=pwr)
 
     def update_pref(self, pref):
-        return pref._replace(p_span0=pref.p0, p_spani=pref.pi - self.loss)
+        return pref._replace(p_span0=pref.p_span0, p_spani=pref.p_spani - self.loss)
 
     def __call__(self, spectral_info):
         carriers = tuple(self.propagate(*spectral_info.carriers))
         pref = self.update_pref(spectral_info.pref)
-        return spectral_info.update(carriers=carriers, pref=pref)
+        return spectral_info._replace(carriers=carriers, pref=pref)
 
 FiberParams = namedtuple('FiberParams', 'type_variety length loss_coef length_units \
                                          att_in con_in con_out dispersion gamma')
@@ -283,12 +291,12 @@ class Fiber(Node):
 
     @property
     def fiber_loss(self):
-        # dB fiber loss, not including padding attenuator
+        """Fiber loss in dB, not including padding attenuator"""
         return self.loss_coef * self.length + self.con_in + self.con_out
 
     @property
     def loss(self):
-        #total loss incluiding padding att_in: useful for polymorphism with roadm loss
+        """total loss including padding att_in: useful for polymorphism with roadm loss"""
         return self.loss_coef * self.length + self.con_in + self.con_out + self.att_in
 
     @property
@@ -313,16 +321,13 @@ class Fiber(Node):
 
     def carriers(self, loc, attr):
         """retrieve carriers information
-        loc = (in, out) of the class element
+
-        attr = (ase, nli, signal, total) power information"""
+        :param loc: (in, out) of the class element
+        :param attr: (ase, nli, signal, total) power information
+        """
         if not (loc in ('in', 'out') and attr in ('nli', 'signal', 'total', 'ase')):
             yield None
             return
-        power_dict = {
-                        'nli':      'nonlinear_interference',
-                        'ase':      'amplified_spontaneous_emission'
-                    }
-        attr = power_dict.get(attr, attr)
         loc_attr = 'carriers_'+loc
         for c in getattr(self, loc_attr) :
             if attr == 'total':
@@ -330,46 +335,30 @@ class Fiber(Node):
             else:
                 yield c.power._asdict().get(attr, None)
 
-    def beta2(self, ref_wavelength=None):
+    def beta2(self, ref_wavelength=1550e-9):
-        """ Returns beta2 from dispersion parameter.
+        """Returns beta2 from dispersion parameter.
         Dispersion is entered in ps/nm/km.
-        Disperion can be a numpy array or a single value.  If a
+        Disperion can be a numpy array or a single value.
-        value ref_wavelength is not entered 1550e-9m will be assumed.
+
-        ref_wavelength can be a numpy array.
+        :param ref_wavelength: can be a numpy array; default: 1550nm
         """
         # TODO|jla: discuss beta2 as method or attribute
-        wl = 1550e-9 if ref_wavelength is None else ref_wavelength
         D = abs(self.dispersion)
-        b2 = (wl ** 2) * D / (2 * pi * c)  # 10^21 scales [ps^2/km]
+        b2 = (ref_wavelength ** 2) * D / (2 * pi * c)  # 10^21 scales [ps^2/km]
         return b2 # s/Hz/m
 
     def dbkm_2_lin(self):
-        """ calculates the linear loss coefficient
+        """calculates the linear loss coefficient"""
-        """
+        # linear loss coefficient in dB/km^-1
-        # alpha_pcoef is linear loss coefficient in dB/km^-1
-        # alpha_acoef is linear loss field amplitude coefficient in m^-1
         alpha_pcoef = self.loss_coef
+        # linear loss field amplitude coefficient in m^-1
         alpha_acoef = alpha_pcoef / (2 * 10 * log10(exp(1)))
         return alpha_pcoef, alpha_acoef
 
-    def _psi(self, carrier, interfering_carrier):
-        """ Calculates eq. 123 from	arXiv:1209.0394.
-        """
-        if carrier.num_chan == interfering_carrier.num_chan: # SCI
-            psi = arcsinh(0.5 * pi**2 * self.asymptotic_length
-                              * abs(self.beta2()) * carrier.baud_rate**2)
-        else: # XCI
-            delta_f = carrier.freq - interfering_carrier.freq
-            psi = arcsinh(pi**2 * self.asymptotic_length * abs(self.beta2())
-                                * carrier.baud_rate * (delta_f + 0.5 * interfering_carrier.baud_rate))
-            psi -= arcsinh(pi**2 * self.asymptotic_length * abs(self.beta2())
-                                 * carrier.baud_rate * (delta_f - 0.5 * interfering_carrier.baud_rate))
-
-        return psi
-
     def _gn_analytic(self, carrier, *carriers):
-        """ Computes the nonlinear interference power on a single carrier.
+        """Computes the nonlinear interference power on a single carrier.
-        The method uses eq. 120 from arXiv:1209.0394.
+        The method uses eq. 120 from `arXiv:1209.0394 <https://arxiv.org/abs/1209.0394>`__.
+
         :param carrier: the signal under analysis
         :param carriers: the full WDM comb
         :return: carrier_nli: the amount of nonlinear interference in W on the under analysis
@@ -377,7 +366,7 @@ class Fiber(Node):
 
         g_nli = 0
         for interfering_carrier in carriers:
-            psi = self._psi(carrier, interfering_carrier)
+            psi = _psi(carrier, interfering_carrier, beta2=self.beta2(), asymptotic_length=self.asymptotic_length)
             g_nli += (interfering_carrier.power.signal/interfering_carrier.baud_rate)**2 \
                      * (carrier.power.signal/carrier.baud_rate) * psi
 
@@ -396,8 +385,8 @@ class Fiber(Node):
         for carrier in carriers:
             pwr = carrier.power
             pwr = pwr._replace(signal=pwr.signal/attenuation,
-                               nonlinear_interference=pwr.nli/attenuation,
+                               nli=pwr.nli/attenuation,
-                               amplified_spontaneous_emission=pwr.ase/attenuation)
+                               ase=pwr.ase/attenuation)
             carrier = carrier._replace(power=pwr)
             chan.append(carrier)
 
@@ -409,20 +398,77 @@ class Fiber(Node):
             pwr = carrier.power
             carrier_nli = self._gn_analytic(carrier, *carriers)
             pwr = pwr._replace(signal=pwr.signal/self.lin_attenuation/attenuation,
-                               nonlinear_interference=(pwr.nli+carrier_nli)/self.lin_attenuation/attenuation,
+                               nli=(pwr.nli+carrier_nli)/self.lin_attenuation/attenuation,
-                               amplified_spontaneous_emission=pwr.ase/self.lin_attenuation/attenuation)
+                               ase=pwr.ase/self.lin_attenuation/attenuation)
             yield carrier._replace(power=pwr)
 
     def update_pref(self, pref):
-        self.pch_out_db = round(pref.pi - self.loss, 2)
+        self.pch_out_db = round(pref.p_spani - self.loss, 2)
-        return pref._replace(p_span0=pref.p0, p_spani=self.pch_out_db)
+        return pref._replace(p_span0=pref.p_span0, p_spani=self.pch_out_db)
 
     def __call__(self, spectral_info):
         self.carriers_in = spectral_info.carriers
         carriers = tuple(self.propagate(*spectral_info.carriers))
         pref = self.update_pref(spectral_info.pref)
         self.carriers_out = carriers
-        return spectral_info.update(carriers=carriers, pref=pref)
+        return spectral_info._replace(carriers=carriers, pref=pref)
+
+RamanFiberParams = namedtuple('RamanFiberParams', 'type_variety length loss_coef length_units \
+                                         att_in con_in con_out dispersion gamma raman_efficiency')
+
+class RamanFiber(Fiber):
+    def __init__(self, *args, params=None, **kwargs):
+        if params is None:
+            params = {}
+        if 'con_in' not in params:
+            # if not defined in the network json connector loss in/out
+            # the None value will be updated in network.py[build_network]
+            # with default values from eqpt_config.json[Spans]
+            params['con_in'] = None
+            params['con_out'] = None
+        if 'att_in' not in params:
+            #fixed attenuator for padding
+            params['att_in'] = 0
+
+        # TODO: can we re-use the Fiber constructor in a better way?
+        Node.__init__(self, *args, params=RamanFiberParams(**params), **kwargs)
+        self.type_variety = self.params.type_variety
+        self.length = self.params.length * UNITS[self.params.length_units] # in m
+        self.loss_coef = self.params.loss_coef * 1e-3 # lineic loss dB/m
+        self.lin_loss_coef = self.params.loss_coef / (20 * log10(exp(1)))
+        self.att_in = self.params.att_in
+        self.con_in = self.params.con_in
+        self.con_out = self.params.con_out
+        self.dispersion = self.params.dispersion  # s/m/m
+        self.gamma = self.params.gamma # 1/W/m
+        self.pch_out_db = None
+        self.carriers_in = None
+        self.carriers_out = None
+        # TODO|jla: discuss factor 2 in the linear lineic attenuation
+
+    @property
+    def sim_params(self):
+        return self._sim_params
+
+    @sim_params.setter
+    def sim_params(self, sim_params=None):
+        self._sim_params = sim_params
+
+    def update_pref(self, pref, *carriers):
+        pch_out_db = lin2db(mean([carrier.power.signal for carrier in carriers])) + 30
+        self.pch_out_db = round(pch_out_db, 2)
+        return pref._replace(p_span0=pref.p_span0, p_spani=self.pch_out_db)
+
+    def __call__(self, spectral_info):
+        self.carriers_in = spectral_info.carriers
+        carriers = tuple(self.propagate(*spectral_info.carriers))
+        pref = self.update_pref(spectral_info.pref, *carriers)
+        self.carriers_out = carriers
+        return spectral_info._replace(carriers=carriers, pref=pref)
+
+    def propagate(self, *carriers):
+        for propagated_carrier in propagate_raman_fiber(self, *carriers):
+            yield propagated_carrier
 
 class EdfaParams:
     def __init__(self, **params):
@@ -469,12 +515,11 @@ class EdfaOperational:
                 f'tilt_target={self.tilt_target!r})')
 
 class Edfa(Node):
-    def __init__(self, *args, params={}, operational={}, **kwargs):
+    def __init__(self, *args, params=None, operational=None, **kwargs):
-        #TBC is this useful? put in comment for now:
+        if params is None:
-        #if params is None:
+            params = {}
-        #    params = {}
+        if operational is None:
-        #if operational is None:
+            operational = {}
-        #    operational = {}
         super().__init__(
             *args,
             params=EdfaParams(**params),
@@ -550,16 +595,13 @@ class Edfa(Node):
 
     def carriers(self, loc, attr):
         """retrieve carriers information
-        loc = (in, out) of the class element
+
-        attr = (ase, nli, signal, total) power information"""
+        :param loc: (in, out) of the class element
+        :param attr: (ase, nli, signal, total) power information
+        """
         if not (loc in ('in', 'out') and attr in ('nli', 'signal', 'total', 'ase')):
             yield None
             return
-        power_dict = {
-                        'nli':      'nonlinear_interference',
-                        'ase':      'amplified_spontaneous_emission'
-                    }
-        attr = power_dict.get(attr, attr)
         loc_attr = 'carriers_'+loc
         for c in getattr(self, loc_attr) :
             if attr == 'total':
@@ -568,12 +610,12 @@ class Edfa(Node):
                 yield c.power._asdict().get(attr, None)
 
     def interpol_params(self, frequencies, pin, baud_rates, pref):
-        """interpolate SI channel frequencies with the edfa dgt and gain_ripple frquencies from json
+        """interpolate SI channel frequencies with the edfa dgt and gain_ripple frquencies from JSON
         set the edfa class __init__ None parameters :
                 self.channel_freq, self.nf, self.interpol_dgt and self.interpol_gain_ripple
         """
         # TODO|jla: read amplifier actual frequencies from additional params in json
-        amplifier_freq = itufl(len(self.params.dgt), self.params.f_min, self.params.f_max) * 1e12 # Hz
+        amplifier_freq = itufl(len(self.params.dgt), self.params.f_min, self.params.f_max) # Hz
         self.channel_freq = frequencies
         self.interpol_dgt = interp(self.channel_freq, amplifier_freq, self.params.dgt)
 
@@ -586,19 +628,19 @@ class Edfa(Node):
         """in power mode: delta_p is defined and can be used to calculate the power target
         This power target is used calculate the amplifier gain"""
         if self.delta_p is not None:
-            self.target_pch_out_db = round(self.delta_p + pref.p0, 2)
+            self.target_pch_out_db = round(self.delta_p + pref.p_span0, 2)
-            self.effective_gain = self.target_pch_out_db - pref.pi
+            self.effective_gain = self.target_pch_out_db - pref.p_spani
 
         """check power saturation and correct effective gain & power accordingly:"""            
         self.effective_gain = min(  
                                     self.effective_gain, 
-                                    self.params.p_max - (pref.pi + pref.neq_ch)
+                                    self.params.p_max - (pref.p_spani + pref.neq_ch)
                                     )
         #print(self.uid, self.effective_gain, self.operational.gain_target)
-        self.effective_pch_out_db = round(pref.pi + self.effective_gain, 2)
+        self.effective_pch_out_db = round(pref.p_spani + self.effective_gain, 2)
 
         """check power saturation and correct target_gain accordingly:"""
-        #print(self.uid, self.effective_gain, self.pin_db, pref.pi)
+        #print(self.uid, self.effective_gain, self.pin_db, pref.p_spani)
         self.nf = self._calc_nf()
         self.gprofile = self._gain_profile(pin)
 
@@ -624,14 +666,8 @@ class Edfa(Node):
             nf_avg = pin_ch - polyval(nf_model.nf_coef, pin_ch) + 58
         elif type_def == 'advanced_model':
             nf_avg = polyval(nf_fit_coeff, -dg)
-        else :
+        else:
-            print(
+            assert False, "Unrecognized amplifier type, this should have been checked by the JSON loader"
-                f'\x1b[1;31;40m'\
-                + f'CRITICAL: unrecognized type def _{self.params.type_def}_\n\
-                    => please check eqpt_config.json'\
-                + '\x1b[0m'
-                )                        
-            exit()            
         return nf_avg+pad, pad
 
     def _calc_nf(self, avg = False):
@@ -673,7 +709,7 @@ class Edfa(Node):
             return self.interpol_nf_ripple + nf_avg # input VOA = 1 for 1 NF degradation
 
     def noise_profile(self, df):
-        """ noise_profile(bw) computes amplifier ase (W) in signal bw (Hz)
+        """noise_profile(bw) computes amplifier ase (W) in signal bw (Hz)
         noise is calculated at amplifier input
 
         :bw: signal bandwidth = baud rate in Hz
@@ -828,7 +864,7 @@ class Edfa(Node):
         return g1st - voa + array(self.interpol_dgt) * dgts3
 
     def propagate(self, pref, *carriers):
-        """add ase noise to the propagating carriers of SpectralInformation"""
+        """add ASE noise to the propagating carriers of :class:`.info.SpectralInformation`"""
         pin = array([c.power.signal+c.power.nli+c.power.ase for c in carriers]) # pin in W
         freq = array([c.frequency for c in carriers])
         brate = array([c.baud_rate for c in carriers])
@@ -842,17 +878,17 @@ class Edfa(Node):
         for gain, carrier_ase, carrier in zip(gains, carrier_ases, carriers):
             pwr = carrier.power
             pwr = pwr._replace(signal=pwr.signal*gain/att,
-                               nonlinear_interference=pwr.nli*gain/att,
+                               nli=pwr.nli*gain/att,
-                               amplified_spontaneous_emission=(pwr.ase+carrier_ase)*gain/att)
+                               ase=(pwr.ase+carrier_ase)*gain/att)
             yield carrier._replace(power=pwr)
 
     def update_pref(self, pref):
-        return pref._replace(p_span0=pref.p0,
+        return pref._replace(p_span0=pref.p_span0,
-                            p_spani=pref.pi + self.effective_gain - self.out_voa)
+                            p_spani=pref.p_spani + self.effective_gain - self.out_voa)
 
     def __call__(self, spectral_info):
         self.carriers_in = spectral_info.carriers
         carriers = tuple(self.propagate(spectral_info.pref, *spectral_info.carriers))
         pref = self.update_pref(spectral_info.pref)
         self.carriers_out = carriers
-        return spectral_info.update(carriers=carriers, pref=pref)
+        return spectral_info._replace(carriers=carriers, pref=pref)

gnpy/core/equipment.py

@@ -9,7 +9,6 @@ This module contains functionality for specifying equipment.
 '''
 
 from numpy import clip, polyval
-from sys import exit
 from operator import itemgetter
 from math import isclose
 from pathlib import Path
@@ -17,6 +16,7 @@ from json import load
 from gnpy.core.utils import lin2db, db2lin, load_json
 from collections import namedtuple
 from gnpy.core.elements import Edfa
+from gnpy.core.exceptions import EquipmentConfigError
 import time
 
 Model_vg = namedtuple('Model_vg', 'nf1 nf2 delta_p')
@@ -27,14 +27,14 @@ Model_dual_stage = namedtuple('Model_dual_stage', 'preamp_variety booster_variet
 
 class common:
     def update_attr(self, default_values, kwargs, name):
-        clean_kwargs = {k:v for k,v in kwargs.items() if v !=''}
+        clean_kwargs = {k:v for k, v in kwargs.items() if v != ''}
-        for k,v in default_values.items():
+        for k, v in default_values.items():
-            setattr(self, k, clean_kwargs.get(k,v))
+            setattr(self, k, clean_kwargs.get(k, v))
-            if k not in clean_kwargs and name != 'Amp' :
+            if k not in clean_kwargs and name != 'Amp':
                 print(f'\x1b[1;31;40m'+
-                    f'\n WARNING missing {k} attribute in eqpt_config.json[{name}]'
+                      f'\n WARNING missing {k} attribute in eqpt_config.json[{name}]'+
-                    f'\n default value is {k} = {v}'
+                      f'\n default value is {k} = {v}'+
-                    + '\x1b[0m')
+                      f'\x1b[0m')
                 time.sleep(1)
 
 class SI(common):
@@ -42,13 +42,13 @@ class SI(common):
     {
         "f_min":            191.35e12,
         "f_max":            196.1e12,
-        "baud_rate":        32e9,        
+        "baud_rate":        32e9,
         "spacing":          50e9,
         "power_dbm":        0,
-        "power_range_db":   [0,0,0.5],
+        "power_range_db":   [0, 0, 0.5],
         "roll_off":         0.15,
         "tx_osnr":          45,
-        "sys_margins":      0    
+        "sys_margins":      0
     }
 
     def __init__(self, **kwargs):
@@ -69,7 +69,7 @@ class Span(common):
         'con_in':                           0,
         'con_out':                          0
     }
-    
+
     def __init__(self, **kwargs):
         self.update_attr(self.default_values, kwargs, 'Span')
 
@@ -77,8 +77,12 @@ class Roadm(common):
     default_values = \
     {
         'target_pch_out_db':   -17,
-        'add_drop_osnr':       100
+        'add_drop_osnr':       100,
-    }    
+        'restrictions': {
+            'preamp_variety_list':[],
+            'booster_variety_list':[]
+            }
+    }
 
     def __init__(self, **kwargs):
         self.update_attr(self.default_values, kwargs, 'Roadm')
@@ -105,6 +109,23 @@ class Fiber(common):
     def __init__(self, **kwargs):
         self.update_attr(self.default_values, kwargs, 'Fiber')
 
+class RamanFiber(common):
+    default_values = \
+    {
+        'type_variety':  '',
+        'dispersion':    None,
+        'gamma':         0,
+        'raman_efficiency': None
+    }
+
+    def __init__(self, **kwargs):
+        self.update_attr(self.default_values, kwargs, 'RamanFiber')
+        for param in ('cr', 'frequency_offset'):
+            if param not in self.raman_efficiency:
+                raise EquipmentConfigError(f'RamanFiber.raman_efficiency: missing "{param}" parameter')
+        if self.raman_efficiency['frequency_offset'] != sorted(self.raman_efficiency['frequency_offset']):
+            raise EquipmentConfigError(f'RamanFiber.raman_efficiency.frequency_offset is not sorted')
+
 class Amp(common):
     default_values = \
     {
@@ -134,7 +155,7 @@ class Amp(common):
         config = Path(filename).parent / 'default_edfa_config.json'
 
         type_variety = kwargs['type_variety']
-        type_def = kwargs.get('type_def', 'variable_gain') #default compatibility with older json eqpt files
+        type_def = kwargs.get('type_def', 'variable_gain') # default compatibility with older json eqpt files
         nf_def = None
         dual_stage_def = None
 
@@ -142,12 +163,12 @@ class Amp(common):
             try:
                 nf0 = kwargs.pop('nf0')
             except KeyError: #nf0 is expected for a fixed gain amp
-                print(f'missing nf0 value input for amplifier: {type_variety} in eqpt_config.json')
+                raise EquipmentConfigError(f'missing nf0 value input for amplifier: {type_variety} in equipment config')
-                exit()
+            for k in ('nf_min', 'nf_max'):
-            try: #remove all remaining nf inputs
+                try:
-                del kwargs['nf_min']
+                    del kwargs[k]
-                del kwargs['nf_max']
+                except KeyError:
-            except KeyError: pass #nf_min and nf_max are not needed for fixed gain amp
+                    pass
             nf_def = Model_fg(nf0)
         elif type_def == 'advanced_model':
             config = Path(filename).parent / kwargs.pop('advanced_config_from_json')
@@ -157,8 +178,7 @@ class Amp(common):
                 nf_min = kwargs.pop('nf_min')
                 nf_max = kwargs.pop('nf_max')
             except KeyError:
-                print(f'missing nf_min/max value input for amplifier: {type_variety} in eqpt_config.json')
+                raise EquipmentConfigError(f'missing nf_min or nf_max value input for amplifier: {type_variety} in equipment config')
-                exit()
             try: #remove all remaining nf inputs
                 del kwargs['nf0']
             except KeyError: pass #nf0 is not needed for variable gain amp
@@ -168,32 +188,28 @@ class Amp(common):
             try:
                 nf_coef = kwargs.pop('nf_coef')
             except KeyError: #nf_coef is expected for openroadm amp
-                print(f'missing nf_coef input for amplifier: {type_variety} in eqpt_config.json')
+                raise EquipmentConfigError(f'missing nf_coef input for amplifier: {type_variety} in equipment config')
-                exit()
             nf_def = Model_openroadm(nf_coef)
         elif type_def == 'dual_stage':
             try: #nf_ram and gain_ram are expected for a hybrid amp
                 preamp_variety = kwargs.pop('preamp_variety')
                 booster_variety = kwargs.pop('booster_variety')
             except KeyError:
-                print(f'missing preamp/booster variety input for amplifier: {type_variety} in eqpt_config.json')
+                raise EquipmentConfigError(f'missing preamp/booster variety input for amplifier: {type_variety} in equipment config')
-                exit()
             dual_stage_def = Model_dual_stage(preamp_variety, booster_variety)
 
         with open(config, encoding='utf-8') as f:
             json_data = load(f)
 
-        return cls(**{**kwargs, **json_data, 
+        return cls(**{**kwargs, **json_data,
             'nf_model': nf_def, 'dual_stage_model': dual_stage_def})
 
 
 def nf_model(type_variety, gain_min, gain_max, nf_min, nf_max):
     if nf_min < -10:
-        print(f'Invalid nf_min value {nf_min!r} for amplifier {type_variety}')
+        raise EquipmentConfigError(f'Invalid nf_min value {nf_min!r} for amplifier {type_variety}')
-        exit()
     if nf_max < -10:
-        print(f'Invalid nf_max value {nf_max!r} for amplifier {type_variety}')
+        raise EquipmentConfigError(f'Invalid nf_max value {nf_max!r} for amplifier {type_variety}')
-        exit()
 
     # NF estimation model based on nf_min and nf_max
     # delta_p:  max power dB difference between first and second stage coils
@@ -208,8 +224,7 @@ def nf_model(type_variety, gain_min, gain_max, nf_min, nf_max):
     nf1 = lin2db(db2lin(nf_min) - db2lin(nf2)/db2lin(g1a_max))
 
     if nf1 < 4:
-        print(f'First coil value too low {nf1} for amplifier {type_variety}')
+        raise EquipmentConfigError(f'First coil value too low {nf1} for amplifier {type_variety}')
-        exit()
 
     # Check 1 dB < delta_p < 6 dB to ensure nf_min and nf_max values make sense.
     # There shouldn't be high nf differences between the two coils:
@@ -221,20 +236,17 @@ def nf_model(type_variety, gain_min, gain_max, nf_min, nf_max):
         delta_p = gain_max - g1a_max
         g1a_min = gain_min - (gain_max-gain_min) - delta_p
         if not 1 < delta_p < 11:
-            print(f'Computed \N{greek capital letter delta}P invalid \
+            raise EquipmentConfigError(f'Computed \N{greek capital letter delta}P invalid \
                 \n 1st coil vs 2nd coil calculated DeltaP {delta_p:.2f} for \
                 \n amplifier {type_variety} is not valid: revise inputs \
                 \n calculated 1st coil NF = {nf1:.2f}, 2nd coil NF = {nf2:.2f}')
-            exit()
     # Check calculated values for nf1 and nf2
     calc_nf_min = lin2db(db2lin(nf1) + db2lin(nf2)/db2lin(g1a_max))
     if not isclose(nf_min, calc_nf_min, abs_tol=0.01):
-        print(f'nf_min does not match calc_nf_min, {nf_min} vs {calc_nf_min} for amp {type_variety}')
+        raise EquipmentConfigError(f'nf_min does not match calc_nf_min, {nf_min} vs {calc_nf_min} for amp {type_variety}')
-        exit()
     calc_nf_max = lin2db(db2lin(nf1) + db2lin(nf2)/db2lin(g1a_min))
     if not isclose(nf_max, calc_nf_max, abs_tol=0.01):
-        print(f'nf_max does not match calc_nf_max, {nf_max} vs {calc_nf_max} for amp {type_variety}')
+        raise EquipmentConfigError(f'nf_max does not match calc_nf_max, {nf_max} vs {calc_nf_max} for amp {type_variety}')
-        exit()
 
     return nf1, nf2, delta_p
 
@@ -256,7 +268,7 @@ def trx_mode_params(equipment, trx_type_variety='', trx_mode='', error_message=F
     """return the trx and SI parameters from eqpt_config for a given type_variety and mode (ie format)"""
     trx_params = {}
     default_si_data = equipment['SI']['default']
-    
+
     try:
         trxs = equipment['Transceiver']
         #if called from path_requests_run.py, trx_mode is filled with None when not specified by user
@@ -269,20 +281,18 @@ def trx_mode_params(equipment, trx_type_variety='', trx_mode='', error_message=F
             trx_params = {**mode_params}
             # sanity check: spacing baudrate must be smaller than min spacing
             if trx_params['baud_rate'] > trx_params['min_spacing'] :
-                msg = f'Inconsistency in equipment library:\n Transpoder "{trx_type_variety}" mode "{trx_params["format"]}" '+\
+                raise EquipmentConfigError(f'Inconsistency in equipment library:\n Transpoder "{trx_type_variety}" mode "{trx_params["format"]}" '+\
-                    f'has baud rate: {trx_params["baud_rate"]*1e-9} GHz greater than min_spacing {trx_params["min_spacing"]*1e-9}.'
+                    f'has baud rate: {trx_params["baud_rate"]*1e-9} GHz greater than min_spacing {trx_params["min_spacing"]*1e-9}.')
-                print(msg)
-                exit()
         else:
             mode_params = {"format": "undetermined",
-                       "baud_rate": None,
+                           "baud_rate": None,
-                       "OSNR": None,
+                           "OSNR": None,
-                       "bit_rate": None,
+                           "bit_rate": None,
-                       "roll_off": None,
+                           "roll_off": None,
-                       "tx_osnr":None,
+                           "tx_osnr":None,
-                       "min_spacing":None,
+                           "min_spacing":None,
-                       "cost":None}
+                           "cost":None}
-            trx_params = {**mode_params} 
+            trx_params = {**mode_params}
         trx_params['f_min'] = equipment['Transceiver'][trx_type_variety].frequency['min']
         trx_params['f_max'] = equipment['Transceiver'][trx_type_variety].frequency['max']
 
@@ -292,9 +302,7 @@ def trx_mode_params(equipment, trx_type_variety='', trx_mode='', error_message=F
         # print(f'spacing {temp}')
     except StopIteration :
         if error_message:
-            print(f'could not find tsp : {trx_type_variety} with mode: {trx_mode} in eqpt library')
+            raise EquipmentConfigError(f'Computation stoped: could not find tsp : {trx_type_variety} with mode: {trx_mode} in eqpt library')
-            print('Computation stopped.')
-            exit()
         else:
             # default transponder charcteristics
             # mainly used with transmission_main_example.py
@@ -311,7 +319,7 @@ def trx_mode_params(equipment, trx_type_variety='', trx_mode='', error_message=F
             nch = automatic_nch(trx_params['f_min'], trx_params['f_max'], trx_params['spacing'])
             trx_params['nb_channel'] = nch
             print(f'There are {nch} channels propagating')
-                
+
     trx_params['power'] =  db2lin(default_si_data.power_dbm)*1e-3
 
     return trx_params
@@ -319,8 +327,8 @@ def trx_mode_params(equipment, trx_type_variety='', trx_mode='', error_message=F
 def automatic_spacing(baud_rate):
     """return the min possible channel spacing for a given baud rate"""
     # TODO : this should parametrized in a cfg file
-    spacing_list = [(33e9,37.5e9), (38e9,50e9), (50e9,62.5e9), (67e9,75e9), (92e9,100e9)] #list of possible tuples
+    # list of possible tuples [(max_baud_rate, spacing_for_this_baud_rate)]
-                                                #[(max_baud_rate, spacing_for_this_baud_rate)]
+    spacing_list = [(33e9, 37.5e9), (38e9, 50e9), (50e9, 62.5e9), (67e9, 75e9), (92e9, 100e9)]
     return min((s[1] for s in spacing_list if s[0] > baud_rate), default=baud_rate*1.2)
 
 def automatic_nch(f_min, f_max, spacing):
@@ -346,24 +354,28 @@ def update_dual_stage(equipment):
         if edfa.type_def == 'dual_stage':
             edfa_preamp = edfa_dict[edfa.dual_stage_model.preamp_variety]
             edfa_booster = edfa_dict[edfa.dual_stage_model.booster_variety]
-            for k,v in edfa_preamp.__dict__.items():
+            for key, value in edfa_preamp.__dict__.items():
-                attr_k = 'preamp_'+k
+                attr_k = 'preamp_' + key
-                setattr(edfa, attr_k, v)
+                setattr(edfa, attr_k, value)
-            for k,v in edfa_booster.__dict__.items():
+            for key, value in edfa_booster.__dict__.items():
-                attr_k = 'booster_'+k
+                attr_k = 'booster_' + key
-                setattr(edfa, attr_k, v)           
+                setattr(edfa, attr_k, value)
             edfa.p_max = edfa_booster.p_max
             edfa.gain_flatmax = edfa_booster.gain_flatmax + edfa_preamp.gain_flatmax
             if edfa.gain_min < edfa_preamp.gain_min:
-                print(
+                raise EquipmentConfigError(f'Dual stage {edfa.type_variety} min gain is lower than its preamp min gain')
-                    f'\x1b[1;31;40m'\
-                    + f'CRITICAL: dual stage {edfa.type_variety} min gain is lower than its preamp min gain\
-                        => please increase its min gain in eqpt_config.json'\
-                    + '\x1b[0m'
-                    )                        
-                exit()
     return equipment
 
+def roadm_restrictions_sanity_check(equipment):
+    """ verifies that booster and preamp restrictions specified in roadm equipment are listed
+    in the edfa.
+    """
+    restrictions = equipment['Roadm']['default'].restrictions['booster_variety_list'] + \
+        equipment['Roadm']['default'].restrictions['preamp_variety_list']
+    for amp_name in restrictions:
+        if amp_name not in equipment['Edfa']:
+            raise EquipmentConfigError(f'ROADM restriction {amp_name} does not refer to a defined EDFA name')
+
 def equipment_from_json(json_data, filename):
     """build global dictionnary eqpt_library that stores all eqpt characteristics:
     edfa type type_variety, fiber type_variety
@@ -378,11 +390,12 @@ def equipment_from_json(json_data, filename):
         equipment[key] = {}
         typ = globals()[key]
         for entry in entries:
-            subkey = entry.get('type_variety', 'default')           
+            subkey = entry.get('type_variety', 'default')
             if key == 'Edfa':
                 equipment[key][subkey] = Amp.from_json(filename, **entry)
-            else:                
+            else:
                 equipment[key][subkey] = typ(**entry)
     equipment = update_trx_osnr(equipment)
     equipment = update_dual_stage(equipment)
+    roadm_restrictions_sanity_check(equipment)
     return equipment

gnpy/core/exceptions.py

@@ -0,0 +1,19 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+'''
+gnpy.core.exceptions
+====================
+
+Exceptions thrown by other gnpy modules
+'''
+
+
+class ConfigurationError(Exception):
+    '''User-provided configuration contains an error'''
+
+class EquipmentConfigError(ConfigurationError):
+    '''Incomplete or wrong configuration within the equipment library'''
+
+class NetworkTopologyError(ConfigurationError):
+    '''Topology of user-provided network is wrong'''

gnpy/core/info.py

@@ -5,7 +5,7 @@
 gnpy.core.info
 ==============
 
-This module contains classes for modelling SpectralInformation.
+This module contains classes for modelling :class:`SpectralInformation`.
 '''
 
 
@@ -16,50 +16,26 @@ from json import loads
 from gnpy.core.utils import load_json
 from gnpy.core.equipment import automatic_nch, automatic_spacing
 
-class ConvenienceAccess:
+class Power(namedtuple('Power', 'signal nli ase')):
-
-    def __init_subclass__(cls):
-        for abbrev, field in getattr(cls, '_ABBREVS', {}).items():
-            setattr(cls, abbrev, property(lambda self, f=field: getattr(self, f)))
-
-    def update(self, **kwargs):
-        for abbrev, field in getattr(self, '_ABBREVS', {}).items():
-            if abbrev in kwargs:
-                kwargs[field] = kwargs.pop(abbrev)
-        return self._replace(**kwargs)
-
-
-class Power(namedtuple('Power', 'signal nonlinear_interference amplified_spontaneous_emission'), ConvenienceAccess):
     """carriers power in W"""
-    _ABBREVS = {'nli': 'nonlinear_interference',
-                'ase': 'amplified_spontaneous_emission',}
 
 
-class Channel(namedtuple('Channel', 'channel_number frequency baud_rate roll_off power'), ConvenienceAccess):
+class Channel(namedtuple('Channel', 'channel_number frequency baud_rate roll_off power')):
+    pass
 
-    _ABBREVS = {'channel':  'channel_number',
-                'num_chan': 'channel_number',
-                'ffs':      'frequency',
-                'freq':     'frequency',}
 
-class Pref(namedtuple('Pref', 'p_span0, p_spani, neq_ch '), ConvenienceAccess):
+class Pref(namedtuple('Pref', 'p_span0, p_spani, neq_ch ')):
     """noiseless reference power in dBm: 
-    p0: inital target carrier power
+    p_span0: inital target carrier power
-    pi: carrier power after element i
+    p_spani: carrier power after element i
-    neqch: equivalent channel count in dB"""
+    neq_ch: equivalent channel count in dB"""
 
-    _ABBREVS = {'p0' :  'p_span0',
-                'pi' :  'p_spani'}
 
-class SpectralInformation(namedtuple('SpectralInformation', 'pref carriers'), ConvenienceAccess):
+class SpectralInformation(namedtuple('SpectralInformation', 'pref carriers')):
 
     def __new__(cls, pref, carriers):
         return super().__new__(cls, pref, carriers)
 
-def merge_input_spectral_information(*si):
-    """mix channel combs of different baud rates and power"""
-    #TODO
-    pass
 
 def create_input_spectral_information(f_min, f_max, roll_off, baud_rate, power, spacing):
     # pref in dB : convert power lin into power in dB
@@ -86,11 +62,11 @@ if __name__ == '__main__':
     si = SpectralInformation()
     spacing = 0.05 # THz
 
-    si = si.update(carriers=tuple(Channel(f+1, 191.3+spacing*(f+1), 32e9, 0.15, Power(1e-3, f, 1)) for f in range(96)))
+    si = si._replace(carriers=tuple(Channel(f+1, 191.3+spacing*(f+1), 32e9, 0.15, Power(1e-3, f, 1)) for f in range(96)))
 
     print(f'si = {si}')
     print(f'si = {si.carriers[0].power.nli}')
     print(f'si = {si.carriers[20].power.nli}')
-    si2 = si.update(carriers=tuple(c.update(power = c.power.update(nli = c.power.nli * 1e5))
+    si2 = si._replace(carriers=tuple(c._replace(power = c.power._replace(nli = c.power.nli * 1e5))
                               for c in si.carriers))
     print(f'si2 = {si2}')

gnpy/core/network.py

@@ -16,11 +16,13 @@ from logging import getLogger
 from os import path
 from operator import itemgetter, attrgetter
 from gnpy.core import elements
-from gnpy.core.elements import Fiber, Edfa, Transceiver, Roadm, Fused
+from gnpy.core.elements import Fiber, Edfa, Transceiver, Roadm, Fused, RamanFiber
 from gnpy.core.equipment import edfa_nf
+from gnpy.core.exceptions import ConfigurationError, NetworkTopologyError
 from gnpy.core.units import UNITS
-from gnpy.core.utils import load_json, save_json, round2float, db2lin, lin2db
+from gnpy.core.utils import (load_json, save_json, round2float, db2lin,
-from sys import exit
+                            merge_amplifier_restrictions)
+from gnpy.core.science_utils import SimParams
 from collections import namedtuple
 
 logger = getLogger(__name__)
@@ -50,13 +52,14 @@ def network_from_json(json_data, equipment):
     for el_config in json_data['elements']:
         typ = el_config.pop('type')
         variety = el_config.pop('type_variety', 'default')
-        if typ in equipment and variety in equipment[typ]:           
+        if typ in equipment and variety in equipment[typ]:
             extra_params = equipment[typ][variety]
-            el_config.setdefault('params', {}).update(extra_params.__dict__)
+            temp = el_config.setdefault('params', {})
-        elif typ in ['Edfa', 'Fiber']: #catch it now because the code will crash later!
+            temp = merge_amplifier_restrictions(temp, extra_params.__dict__)
-            print( f'The {typ} of variety type {variety} was not recognized:'
+            el_config['params'] = temp
+        elif typ in ['Edfa', 'Fiber']: # catch it now because the code will crash later!
+            raise ConfigurationError(f'The {typ} of variety type {variety} was not recognized:'
                     '\nplease check it is properly defined in the eqpt_config json file')
-            exit()
         cls = getattr(elements, typ)
         el = cls(**el_config)
         g.add_node(el)
@@ -66,17 +69,13 @@ def network_from_json(json_data, equipment):
     for cx in json_data['connections']:
         from_node, to_node = cx['from_node'], cx['to_node']
         try:
-            if isinstance(nodes[from_node], Fiber): 
+            if isinstance(nodes[from_node], Fiber):
                 edge_length = nodes[from_node].params.length
-                # print(from_node)
-                # print(edge_length)
             else:
                 edge_length = 0.01
             g.add_edge(nodes[from_node], nodes[to_node], weight = edge_length)
         except KeyError:
-            msg = f'In {__name__} network_from_json function:\n\tcan not find {from_node} or {to_node} defined in {cx}'
+            raise NetworkTopologyError(f'can not find {from_node} or {to_node} defined in {cx}')
-            print(msg)
-            exit(1)
 
     return g
 
@@ -93,21 +92,27 @@ def network_to_json(network):
     data.update(connections)
     return data
 
-def select_edfa(raman_allowed, gain_target, power_target, equipment, uid):
+def select_edfa(raman_allowed, gain_target, power_target, equipment, uid, restrictions=None):
     """amplifer selection algorithm
     @Orange Jean-Luc Augé
     """
     Edfa_list = namedtuple('Edfa_list', 'variety power gain_min nf')
     TARGET_EXTENDED_GAIN = equipment['Span']['default'].target_extended_gain
-    edfa_dict = equipment['Edfa']
+
+    # for roadm restriction only: create a dict including not allowed for design amps
+    # because main use case is to have specific radm amp which are not allowed for ILA
+    # with the auto design
+    edfa_dict = {name: amp for (name, amp) in equipment['Edfa'].items()
+        if restrictions is None or name in restrictions}
+
     pin = power_target - gain_target
 
-    #create 2 list of available amplifiers with relevant attributs for their selection
+    # create 2 list of available amplifiers with relevant attributes for their selection
 
-    #edfa list with :
+    # edfa list with:
-    #extended gain min allowance of 3dB: could be parametrized, but a bit complex
+    # extended gain min allowance of 3dB: could be parametrized, but a bit complex
-    #extended gain max allowance TARGET_EXTENDED_GAIN is coming from eqpt_config.json
+    # extended gain max allowance TARGET_EXTENDED_GAIN is coming from eqpt_config.json
-    #power attribut include power AND gain limitations
+    # power attribut include power AND gain limitations
     edfa_list = [Edfa_list(
                 variety=edfa_variety,
                 power=min(
@@ -122,7 +127,7 @@ def select_edfa(raman_allowed, gain_target, power_target, equipment, uid):
                     -edfa.gain_min,
                 nf=edfa_nf(gain_target, edfa_variety, equipment)) \
                 for edfa_variety, edfa in edfa_dict.items()
-                if (edfa.allowed_for_design and not edfa.raman)]
+                if ((edfa.allowed_for_design or restrictions is not None) and not edfa.raman)]
 
     #consider a Raman list because of different gain_min requirement: 
     #do not allow extended gain min for Raman
@@ -156,15 +161,11 @@ def select_edfa(raman_allowed, gain_target, power_target, equipment, uid):
         #and raman padding at the amplifier input is impossible!
 
         if len(edfa_list) < 1:
-            print(
+            raise ConfigurationError(f'auto_design could not find any amplifier \
-                f'\x1b[1;31;40m'\
-                + f'CRITICAL _ ABORT: auto_design could not find any amplifier \
                     to satisfy min gain requirement in node {uid} \
-                    please increase span fiber padding'\
+                    please increase span fiber padding')
-                + '\x1b[0m'
-                )
-            exit()
         else:
+            # TODO: convert to logging
             print(
                 f'\x1b[1;31;40m'\
                 + f'WARNING: target gain in node {uid} is below all available amplifiers min gain: \
@@ -216,9 +217,8 @@ def target_power(network, node, equipment): #get_fiber_dp
         dp = max(dp_range[0], dp)
         dp = min(dp_range[1], dp)
     except KeyError:
-        print(f'invalid delta_power_range_db definition in eqpt_config[Span]'
+        raise ConfigurationError(f'invalid delta_power_range_db definition in eqpt_config[Span]'
               f'delta_power_range_db: [lower_bound, upper_bound, step]')
-        exit()
 
     if isinstance(node, Roadm):
         dp = 0
@@ -231,10 +231,7 @@ def prev_node_generator(network, node):
     try:
         prev_node = next(n for n in network.predecessors(node))
     except StopIteration:
-        msg = f'In {__name__} prev_node_generator function:\n\t{node.uid} is not properly connected, please check network topology'
+        raise NetworkTopologyError(f'Node {node.uid} is not properly connected, please check network topology')
-        print(msg)
-        logger.critical(msg)
-        exit(1)
     # yield and re-iterate
     if isinstance(prev_node, Fused) or isinstance(node, Fused):
         yield prev_node
@@ -248,8 +245,7 @@ def next_node_generator(network, node):
     try:
         next_node = next(n for n in network.successors(node))
     except StopIteration:
-        print(f'In {__name__} next_node_generator function:\n\t{node.uid}  is not properly connected, please check network topology')
+        raise NetworkTopologyError('Node {node.uid} is not properly connected, please check network topology')
-        exit(1)        
     # yield and re-iterate
     if isinstance(next_node, Fused) or isinstance(node, Fused):
         yield next_node
@@ -338,7 +334,7 @@ def set_egress_amplifier(network, roadm, equipment, pref_total_db):
                 else: #gain mode with effective_gain 
                     gain_target = node.effective_gain
                     dp = prev_dp - node_loss + gain_target
-                #print(node.delta_p, dp, gain_target)
+
                 power_target = pref_total_db + dp         
 
                 raman_allowed = False
@@ -347,9 +343,24 @@ def set_egress_amplifier(network, roadm, equipment, pref_total_db):
                             equipment['Span']['default'].max_fiber_lineic_loss_for_raman
                     raman_allowed = prev_node.params.loss_coef < max_fiber_lineic_loss_for_raman
 
-                if node.params.type_variety == '' :                   
+                # implementation of restrictions on roadm boosters
+                if isinstance(prev_node,Roadm):
+                    if prev_node.restrictions['booster_variety_list']:
+                        restrictions = prev_node.restrictions['booster_variety_list']
+                    else:
+                        restrictions = None
+                elif isinstance(next_node,Roadm):
+                    # implementation of restrictions on roadm preamp
+                    if next_node.restrictions['preamp_variety_list']:
+                        restrictions = next_node.restrictions['preamp_variety_list']
+                    else:
+                        restrictions = None
+                else:
+                    restrictions = None
+
+                if node.params.type_variety == '':                   
                     edfa_variety, power_reduction = select_edfa(raman_allowed, 
-                                   gain_target, power_target, equipment, node.uid)
+                                   gain_target, power_target, equipment, node.uid, restrictions)
                     extra_params = equipment['Edfa'][edfa_variety]
                     node.params.update_params(extra_params.__dict__)
                     dp += power_reduction
@@ -363,7 +374,7 @@ def set_egress_amplifier(network, roadm, equipment, pref_total_db):
                         )                    
                                 
                 node.delta_p = dp if power_mode else None
-                node.effective_gain = gain_target                    
+                node.effective_gain = gain_target
                 set_amplifier_voa(node, power_target, power_mode)
             if isinstance(next_node, Roadm) or isinstance(next_node, Transceiver):
                 break
@@ -438,9 +449,7 @@ def split_fiber(network, fiber, bounds, target_length, equipment):
         next_node = next(network.successors(fiber))
         prev_node = next(network.predecessors(fiber))
     except StopIteration:
-
+        raise NetworkTopologyError(f'Fiber {fiber.uid} is not properly connected, please check network topology')
-        print(f'In {__name__} split_fiber function:\n\t{fiber.uid}   is not properly connected, please check network topology')
-        exit()
 
     network.remove_node(fiber)
 
@@ -493,18 +502,18 @@ def add_fiber_padding(network, fibers, padding):
         try:
             next_node = next(network.successors(fiber))
         except StopIteration:
-            msg = f'In {__name__} add_fiber_padding function:\n\t{fiber.uid}   is not properly connected, please check network topology'
+            raise NetworkTopologyError(f'Fiber {fiber.uid} is not properly connected, please check network topology')
-            print(msg)
-            logger.critical(msg)
-            exit(1)            
         if this_span_loss < padding and not (isinstance(next_node, Fused)):
             #add a padding att_in at the input of the 1st fiber:
             #address the case when several fibers are spliced together
             first_fiber = find_first_node(network, fiber)
-            if first_fiber.att_in is None:
+            # in order to support no booster , fused might be placed
-                first_fiber.att_in = padding - this_span_loss
+            # just after a roadm: need to check that first_fiber is really a fiber
-            else :
+            if isinstance(first_fiber,Fiber):
-                first_fiber.att_in = first_fiber.att_in + padding - this_span_loss
+                if first_fiber.att_in is None:
+                    first_fiber.att_in = padding - this_span_loss
+                else:
+                    first_fiber.att_in = first_fiber.att_in + padding - this_span_loss
 
 def build_network(network, equipment, pref_ch_db, pref_total_db):
     default_span_data = equipment['Span']['default']
@@ -527,6 +536,7 @@ def build_network(network, equipment, pref_ch_db, pref_total_db):
 
     amplified_nodes = [n for n in network.nodes()
                         if isinstance(n, Fiber) or isinstance(n, Roadm)]
+
     for node in amplified_nodes:
         add_egress_amplifier(network, node)
 
@@ -540,3 +550,11 @@ def build_network(network, equipment, pref_ch_db, pref_total_db):
         for t in trx:
             set_egress_amplifier(network, t, equipment, pref_total_db)
 
+def load_sim_params(filename):
+    sim_params = load_json(filename)
+    return SimParams(params=sim_params)
+
+def configure_network(network, sim_params):
+    for node in network.nodes:
+        if isinstance(node, RamanFiber):
+            node.sim_params = sim_params

gnpy/core/node.py

@@ -8,13 +8,13 @@ gnpy.core.node
 This module contains the base class for a network element.
 
 Strictly, a network element is any callable which accepts an immutable
-.info.SpectralInformation object and returns a .info.SpectralInformation object
+:class:`.info.SpectralInformation` object and returns an :class:`.info.SpectralInformation` object
-(a copy.)
+(a copy).
 
 Network elements MUST implement two attributes .uid and .name representing a
 unique identifier and a printable name.
 
-This base class provides a mode convenient way to define a network element
+This base class provides a more convenient way to define a network element
 via subclassing.
 '''
 
@@ -26,10 +26,12 @@ class Location(namedtuple('Location', 'latitude longitude city region')):
         return super().__new__(cls, latitude, longitude, city, region)
 
 class Node:
-    def __init__(self, uid, name=None, params=None, metadata={'location':{}}, operational=None):
+    def __init__(self, uid, name=None, params=None, metadata=None, operational=None):
         if name is None:
             name = uid
         self.uid, self.name = uid, name
+        if metadata is None:
+            metadata = {'location': {}}
         if metadata and not isinstance(metadata.get('location'), Location):
             metadata['location'] = Location(**metadata.pop('location', {}))
         self.params, self.metadata, self.operational = params, metadata, operational

gnpy/core/request.py

@@ -393,18 +393,16 @@ def compute_constrained_path(network, req):
 
     return total_path
 
-def propagate(path, req, equipment, show=False):
+def propagate(path, req, equipment):
     si = create_input_spectral_information(
         req.f_min, req.f_max, req.roll_off, req.baud_rate,
         req.power, req.spacing)
     for el in path:
         si = el(si)
-        if show :
-            print(el)
     path[-1].update_snr(req.tx_osnr, equipment['Roadm']['default'].add_drop_osnr)
     return path
 
-def propagate2(path, req, equipment, show=False):
+def propagate2(path, req, equipment):
     si = create_input_spectral_information(
         req.f_min, req.f_max, req.roll_off, req.baud_rate,
         req.power, req.spacing)
@@ -413,12 +411,10 @@ def propagate2(path, req, equipment, show=False):
         before_si = si
         after_si  = si = el(si)
         infos[el] = before_si, after_si
-        if show :
-            print(el)
     path[-1].update_snr(req.tx_osnr, equipment['Roadm']['default'].add_drop_osnr)
     return infos
 
-def propagate_and_optimize_mode(path, req, equipment, show=False):
+def propagate_and_optimize_mode(path, req, equipment):
     # if mode is unknown : loops on the modes starting from the highest baudrate fiting in the
     # step 1: create an ordered list of modes based on baudrate
     baudrate_to_explore = list(set([m['baud_rate'] for m in equipment['Transceiver'][req.tsp].mode 
@@ -442,8 +438,6 @@ def propagate_and_optimize_mode(path, req, equipment, show=False):
             b, req.power, req.spacing)
             for el in path:
                 si = el(si)
-                if show:
-                    print(el)
             for m in modes_to_explore :
                 if path[-1].snr is not None:
                     path[-1].update_snr(m['tx_osnr'], equipment['Roadm']['default'].add_drop_osnr)

gnpy/core/science_utils.py

@@ -0,0 +1,820 @@
+import numpy as np
+from operator import attrgetter
+from collections import namedtuple
+from logging import getLogger
+import scipy.constants as ph
+from scipy.integrate import solve_bvp
+from scipy.integrate import cumtrapz
+from scipy.interpolate import interp1d
+from scipy.optimize import OptimizeResult
+
+from gnpy.core.utils import db2lin
+
+
+logger = getLogger(__name__)
+
+
+class RamanParams():
+    def __init__(self, params):
+        self._flag_raman = params['flag_raman']
+        self._space_resolution = params['space_resolution']
+        self._tolerance = params['tolerance']
+
+    @property
+    def flag_raman(self):
+        return self._flag_raman
+
+    @property
+    def space_resolution(self):
+        return self._space_resolution
+
+    @property
+    def tolerance(self):
+        return self._tolerance
+
+class NLIParams():
+    def __init__(self, params):
+        self._nli_method_name = params['nli_method_name']
+        self._wdm_grid_size = params['wdm_grid_size']
+        self._dispersion_tolerance = params['dispersion_tolerance']
+        self._phase_shift_tollerance = params['phase_shift_tollerance']
+        self._f_cut_resolution = None
+        self._f_pump_resolution = None
+
+    @property
+    def nli_method_name(self):
+        return self._nli_method_name
+
+    @property
+    def wdm_grid_size(self):
+        return self._wdm_grid_size
+
+    @property
+    def dispersion_tolerance(self):
+        return self._dispersion_tolerance
+
+    @property
+    def phase_shift_tollerance(self):
+        return self._phase_shift_tollerance
+
+    @property
+    def f_cut_resolution(self):
+        return self._f_cut_resolution
+
+    @f_cut_resolution.setter
+    def f_cut_resolution(self, f_cut_resolution):
+        self._f_cut_resolution = f_cut_resolution
+
+    @property
+    def f_pump_resolution(self):
+        return self._f_pump_resolution
+
+    @f_pump_resolution.setter
+    def f_pump_resolution(self, f_pump_resolution):
+        self._f_pump_resolution = f_pump_resolution
+
+class SimParams():
+    def __init__(self, params):
+        self._raman_computed_channels = params['raman_computed_channels']
+        self._raman_params = RamanParams(params=params['raman_parameters'])
+        self._nli_params = NLIParams(params=params['nli_parameters'])
+
+    @property
+    def raman_computed_channels(self):
+        return self._raman_computed_channels
+
+    @property
+    def raman_params(self):
+        return self._raman_params
+
+    @property
+    def nli_params(self):
+        return self._nli_params
+
+class FiberParams():
+    def __init__(self, fiber):
+        self._loss_coef = 2 * fiber.dbkm_2_lin()[1]
+        self._length = fiber.length
+        self._gamma = fiber.gamma
+        self._beta2 = fiber.beta2()
+        self._beta3 = fiber.beta3 if hasattr(fiber, 'beta3') else 0
+        self._f_ref_beta = fiber.f_ref_beta if hasattr(fiber, 'f_ref_beta') else 0
+        self._raman_efficiency = fiber.params.raman_efficiency
+        self._temperature = fiber.operational['temperature']
+
+    @property
+    def loss_coef(self):
+        return self._loss_coef
+
+    @property
+    def length(self):
+        return self._length
+
+    @property
+    def gamma(self):
+        return self._gamma
+
+    @property
+    def beta2(self):
+        return self._beta2
+
+    @property
+    def beta3(self):
+        return self._beta3
+
+    @property
+    def f_ref_beta(self):
+        return self._f_ref_beta
+
+    @property
+    def raman_efficiency(self):
+        return self._raman_efficiency
+
+    @property
+    def temperature(self):
+        return self._temperature
+
+    def alpha0(self, f_ref=193.5e12):
+        """ It returns the zero element of the series expansion of attenuation coefficient alpha(f) in the
+        reference frequency f_ref
+
+        :param f_ref: reference frequency of series expansion [Hz]
+        :return: alpha0: power attenuation coefficient in f_ref [Neper/m]
+        """
+        if not hasattr(self.loss_coef, 'alpha_power'):
+            alpha0 = self.loss_coef
+        else:
+            alpha_interp = interp1d(self.loss_coef['frequency'],
+                                    self.loss_coef['alpha_power'])
+            alpha0 = alpha_interp(f_ref)
+        return alpha0
+
+pump = namedtuple('RamanPump', 'power frequency propagation_direction')
+
+def propagate_raman_fiber(fiber, *carriers):
+    sim_params = fiber.sim_params
+    raman_params = fiber.sim_params.raman_params
+    nli_params = fiber.sim_params.nli_params
+    # apply input attenuation to carriers
+    attenuation_in = db2lin(fiber.con_in + fiber.att_in)
+    chan = []
+    for carrier in carriers:
+        pwr = carrier.power
+        pwr = pwr._replace(signal=pwr.signal / attenuation_in,
+                           nli=pwr.nli / attenuation_in,
+                           ase=pwr.ase / attenuation_in)
+        carrier = carrier._replace(power=pwr)
+        chan.append(carrier)
+    carriers = tuple(f for f in chan)
+    fiber_params = FiberParams(fiber)
+
+    # evaluate fiber attenuation involving also SRS if required by sim_params
+    if 'raman_pumps' in fiber.operational:
+        raman_pumps = tuple(pump(p['power'], p['frequency'], p['propagation_direction'])
+                            for p in fiber.operational['raman_pumps'])
+    else:
+        raman_pumps = None
+    raman_solver = RamanSolver(raman_params=raman_params, fiber_params=fiber_params)
+    stimulated_raman_scattering = raman_solver.stimulated_raman_scattering(carriers=carriers,
+                                                                           raman_pumps=raman_pumps)
+    fiber_attenuation = (stimulated_raman_scattering.rho[:, -1])**-2
+    if not raman_params.flag_raman:
+        fiber_attenuation = tuple(fiber.lin_attenuation for _ in carriers)
+
+    # evaluate Raman ASE noise if required by sim_params and if raman pumps are present
+    if raman_params.flag_raman and raman_pumps:
+        raman_ase = raman_solver.spontaneous_raman_scattering.power[:, -1]
+    else:
+        raman_ase = tuple(0 for _ in carriers)
+
+    # evaluate nli and propagate in fiber
+    attenuation_out = db2lin(fiber.con_out)
+    nli_solver = NliSolver(nli_params=nli_params, fiber_params=fiber_params)
+    nli_solver.stimulated_raman_scattering = stimulated_raman_scattering
+
+    nli_frequencies = []
+    computed_nli = []
+    for carrier in (c for c in carriers if c.channel_number in sim_params.raman_computed_channels):
+        resolution_param = frequency_resolution(carrier, carriers, sim_params, fiber_params)
+        f_cut_resolution, f_pump_resolution, _, _ = resolution_param
+        nli_params.f_cut_resolution = f_cut_resolution
+        nli_params.f_pump_resolution = f_pump_resolution
+        nli_frequencies.append(carrier.frequency)
+        computed_nli.append(nli_solver.compute_nli(carrier, *carriers))
+
+    new_carriers = []
+    for carrier, attenuation, rmn_ase in zip(carriers, fiber_attenuation, raman_ase):
+        carrier_nli = np.interp(carrier.frequency, nli_frequencies, computed_nli)
+        pwr = carrier.power
+        pwr = pwr._replace(signal=pwr.signal/attenuation/attenuation_out,
+                           nli=(pwr.nli+carrier_nli)/attenuation/attenuation_out,
+                           ase=((pwr.ase/attenuation)+rmn_ase)/attenuation_out)
+        new_carriers.append(carrier._replace(power=pwr))
+    return new_carriers
+
+def frequency_resolution(carrier, carriers, sim_params, fiber_params):
+    def _get_freq_res_k_phi(delta_count, grid_size, alpha0, delta_z, beta2, k_tol, phi_tol):
+        res_phi = _get_freq_res_phase_rotation(delta_count, grid_size, delta_z, beta2, phi_tol)
+        res_k = _get_freq_res_dispersion_attenuation(delta_count, grid_size, alpha0, beta2, k_tol)
+        res_dict = {'res_phi': res_phi, 'res_k': res_k}
+        method = min(res_dict, key=res_dict.get)
+        return res_dict[method], method, res_dict
+
+    def _get_freq_res_dispersion_attenuation(delta_count, grid_size, alpha0, beta2, k_tol):
+        return k_tol * abs(alpha0) / abs(beta2) / (1 + delta_count) / (4 * np.pi ** 2 * grid_size)
+
+    def _get_freq_res_phase_rotation(delta_count, grid_size, delta_z, beta2, phi_tol):
+        return phi_tol / abs(beta2) / (1 + delta_count) / delta_z / (4 * np.pi ** 2 * grid_size)
+
+    grid_size = sim_params.nli_params.wdm_grid_size
+    delta_z = sim_params.raman_params.space_resolution
+    alpha0 = fiber_params.alpha0()
+    beta2 = fiber_params.beta2
+    k_tol = sim_params.nli_params.dispersion_tolerance
+    phi_tol = sim_params.nli_params.phase_shift_tollerance
+    f_pump_resolution, method_f_pump, res_dict_pump = \
+        _get_freq_res_k_phi(0, grid_size, alpha0, delta_z, beta2, k_tol, phi_tol)
+    f_cut_resolution = {}
+    method_f_cut = {}
+    res_dict_cut = {}
+    for cut_carrier in carriers:
+        delta_number = cut_carrier.channel_number - carrier.channel_number
+        delta_count = abs(delta_number)
+        f_res, method, res_dict = \
+            _get_freq_res_k_phi(delta_count, grid_size, alpha0, delta_z, beta2, k_tol, phi_tol)
+        f_cut_resolution[f'delta_{delta_number}'] = f_res
+        method_f_cut[delta_number] = method
+        res_dict_cut[delta_number] = res_dict
+    return [f_cut_resolution, f_pump_resolution, (method_f_cut, method_f_pump), (res_dict_cut, res_dict_pump)]
+
+def raised_cosine_comb(f, *carriers):
+    """ Returns an array storing the PSD of a WDM comb of raised cosine shaped
+    channels at the input frequencies defined in array f
+    :param f: numpy array of frequencies in Hz
+    :param carriers: namedtuple describing the WDM comb
+    :return: PSD of the WDM comb evaluated over f
+    """
+    psd = np.zeros(np.shape(f))
+    for carrier in carriers:
+        f_nch = carrier.frequency
+        g_ch = carrier.power.signal / carrier.baud_rate
+        ts = 1 / carrier.baud_rate
+        passband = (1 - carrier.roll_off) / (2 / carrier.baud_rate)
+        stopband = (1 + carrier.roll_off) / (2 / carrier.baud_rate)
+        ff = np.abs(f - f_nch)
+        tf = ff - passband
+        if carrier.roll_off == 0:
+            psd = np.where(tf <= 0, g_ch, 0.) + psd
+        else:
+            psd = g_ch * (np.where(tf <= 0, 1., 0.) + 1 / 2 * (1 + np.cos(np.pi * ts / carrier.roll_off * tf)) *
+                          np.where(tf > 0, 1., 0.) * np.where(np.abs(ff) <= stopband, 1., 0.)) + psd
+    return psd
+
+class RamanSolver:
+    def __init__(self, raman_params=None, fiber_params=None):
+        """ Initialize the fiber object with its physical parameters
+        :param length: fiber length in m.
+        :param alphap: fiber power attenuation coefficient vs frequency in 1/m. numpy array
+        :param freq_alpha: frequency axis of alphap in Hz. numpy array
+        :param cr_raman: Raman efficiency vs frequency offset in 1/W/m. numpy array
+        :param freq_cr: reference frequency offset axis for cr_raman. numpy array
+        :param raman_params: namedtuple containing the solver parameters (optional).
+        """
+        self.fiber_params = fiber_params
+        self.raman_params = raman_params
+        self._carriers = None
+        self._stimulated_raman_scattering = None
+        self._spontaneous_raman_scattering = None
+
+    @property
+    def fiber_params(self):
+        return self._fiber_params
+
+    @fiber_params.setter
+    def fiber_params(self, fiber_params):
+        self._stimulated_raman_scattering = None
+        self._fiber_params = fiber_params
+
+    @property
+    def carriers(self):
+        return self._carriers
+
+    @carriers.setter
+    def carriers(self, carriers):
+        """
+        :param carriers: tuple of namedtuples containing information about carriers
+        :return:
+        """
+        self._carriers = carriers
+        self._stimulated_raman_scattering = None
+
+    @property
+    def raman_pumps(self):
+        return self._raman_pumps
+
+    @raman_pumps.setter
+    def raman_pumps(self, raman_pumps):
+        self._raman_pumps = raman_pumps
+        self._stimulated_raman_scattering = None
+
+    @property
+    def raman_params(self):
+        return self._raman_params
+
+    @raman_params.setter
+    def raman_params(self, raman_params):
+        """
+        :param raman_params: namedtuple containing the solver parameters (optional).
+        :return:
+        """
+        self._raman_params = raman_params
+        self._stimulated_raman_scattering = None
+        self._spontaneous_raman_scattering = None
+
+    @property
+    def spontaneous_raman_scattering(self):
+        if self._spontaneous_raman_scattering is None:
+            # SET STUFF
+            loss_coef = self.fiber_params.loss_coef
+            raman_efficiency = self.fiber_params.raman_efficiency
+            temperature = self.fiber_params.temperature
+            carriers = self.carriers
+            raman_pumps = self.raman_pumps
+
+            logger.debug('Start computing fiber Spontaneous Raman Scattering')
+            power_spectrum, freq_array, prop_direct, bn_array = self._compute_power_spectrum(carriers, raman_pumps)
+
+            if not hasattr(loss_coef, 'alpha_power'):
+                alphap_fiber = loss_coef * np.ones(freq_array.shape)
+            else:
+                interp_alphap = interp1d(loss_coef['frequency'], loss_coef['alpha_power'])
+                alphap_fiber = interp_alphap(freq_array)
+
+            freq_diff = abs(freq_array - np.reshape(freq_array, (len(freq_array), 1)))
+            interp_cr = interp1d(raman_efficiency['frequency_offset'], raman_efficiency['cr'])
+            cr = interp_cr(freq_diff)
+
+            # z propagation axis
+            z_array = self._stimulated_raman_scattering.z
+            ase_bc = np.zeros(freq_array.shape)
+
+            # calculate ase power
+            spontaneous_raman_scattering = self._int_spontaneous_raman(z_array, self._stimulated_raman_scattering.power,
+                                                                       alphap_fiber, freq_array, cr, freq_diff, ase_bc,
+                                                                       bn_array, temperature)
+
+            setattr(spontaneous_raman_scattering, 'frequency', freq_array)
+            setattr(spontaneous_raman_scattering, 'z', z_array)
+            setattr(spontaneous_raman_scattering, 'power', spontaneous_raman_scattering.x)
+            delattr(spontaneous_raman_scattering, 'x')
+
+            logger.debug(spontaneous_raman_scattering.message)
+
+            self._spontaneous_raman_scattering = spontaneous_raman_scattering
+
+        return self._spontaneous_raman_scattering
+
+    @staticmethod
+    def _compute_power_spectrum(carriers, raman_pumps=None):
+        """
+        Rearrangement of spectral and Raman pump information to make them compatible with Raman solver
+        :param carriers: a tuple of namedtuples describing the transmitted channels
+        :param raman_pumps: a namedtuple describing the Raman pumps
+        :return:
+        """
+
+        # Signal power spectrum
+        pow_array = np.array([])
+        f_array = np.array([])
+        noise_bandwidth_array = np.array([])
+        for carrier in sorted(carriers, key=attrgetter('frequency')):
+            f_array = np.append(f_array, carrier.frequency)
+            pow_array = np.append(pow_array, carrier.power.signal)
+            ref_bw = carrier.baud_rate
+            noise_bandwidth_array = np.append(noise_bandwidth_array, ref_bw)
+
+        propagation_direction = np.ones(len(f_array))
+
+        # Raman pump power spectrum
+        if raman_pumps:
+            for pump in raman_pumps:
+                pow_array = np.append(pow_array, pump.power)
+                f_array = np.append(f_array, pump.frequency)
+                direction = +1 if pump.propagation_direction.lower() == 'coprop' else -1
+                propagation_direction = np.append(propagation_direction, direction)
+                noise_bandwidth_array = np.append(noise_bandwidth_array, ref_bw)
+
+        # Final sorting
+        ind = np.argsort(f_array)
+        f_array = f_array[ind]
+        pow_array = pow_array[ind]
+        propagation_direction = propagation_direction[ind]
+
+        return pow_array, f_array, propagation_direction, noise_bandwidth_array
+
+    def _int_spontaneous_raman(self, z_array, raman_matrix, alphap_fiber, freq_array, cr_raman_matrix, freq_diff, ase_bc, bn_array, temperature):
+        spontaneous_raman_scattering = OptimizeResult()
+
+        dx = self.raman_params.space_resolution
+        h = ph.value('Planck constant')
+        kb = ph.value('Boltzmann constant')
+
+        power_ase = np.nan * np.ones(raman_matrix.shape)
+        int_pump = cumtrapz(raman_matrix, z_array, dx=dx, axis=1, initial=0)
+
+        for f_ind, f_ase in enumerate(freq_array):
+            cr_raman = cr_raman_matrix[f_ind, :]
+            vibrational_loss = f_ase / freq_array[:f_ind]
+            eta = 1/(np.exp((h*freq_diff[f_ind, f_ind+1:])/(kb*temperature)) - 1)
+
+            int_fiber_loss = -alphap_fiber[f_ind] * z_array
+            int_raman_loss = np.sum((cr_raman[:f_ind] * vibrational_loss * int_pump[:f_ind, :].transpose()).transpose(), axis=0)
+            int_raman_gain = np.sum((cr_raman[f_ind + 1:] * int_pump[f_ind + 1:, :].transpose()).transpose(), axis=0)
+
+            int_gain_loss = int_fiber_loss + int_raman_gain + int_raman_loss
+
+            new_ase = np.sum((cr_raman[f_ind+1:] * (1 + eta) * raman_matrix[f_ind+1:, :].transpose()).transpose() * h * f_ase * bn_array[f_ind], axis=0)
+
+            bc_evolution = ase_bc[f_ind] * np.exp(int_gain_loss)
+            ase_evolution = np.exp(int_gain_loss) * cumtrapz(new_ase*np.exp(-int_gain_loss), z_array, dx=dx, initial=0)
+
+            power_ase[f_ind, :] = bc_evolution + ase_evolution
+
+        spontaneous_raman_scattering.x = 2 * power_ase
+        spontaneous_raman_scattering.success = True
+        spontaneous_raman_scattering.message = "Spontaneous Raman Scattering evaluated successfully"
+
+        return spontaneous_raman_scattering
+
+    def stimulated_raman_scattering(self, carriers, raman_pumps=None):
+        """ Returns stimulated Raman scattering solution including 
+        fiber gain/loss profile.
+        :return: self._stimulated_raman_scattering: the SRS problem solution.
+        scipy.interpolate.PPoly instance
+        """
+
+        if self._stimulated_raman_scattering is None:
+            # fiber parameters
+            fiber_length = self.fiber_params.length
+            loss_coef = self.fiber_params.loss_coef
+            if self.raman_params.flag_raman:
+                raman_efficiency = self.fiber_params.raman_efficiency
+            else:
+                raman_efficiency = self.fiber_params.raman_efficiency
+                raman_efficiency['cr'] = np.array(raman_efficiency['cr']) * 0
+            # raman solver parameters
+            z_resolution = self.raman_params.space_resolution
+            tolerance = self.raman_params.tolerance
+
+            logger.debug('Start computing fiber Stimulated Raman Scattering')
+
+            power_spectrum, freq_array, prop_direct, _ = self._compute_power_spectrum(carriers, raman_pumps)
+
+            if not hasattr(loss_coef, 'alpha_power'):
+                alphap_fiber = loss_coef * np.ones(freq_array.shape)
+            else:
+                interp_alphap = interp1d(loss_coef['frequency'], loss_coef['alpha_power'])
+                alphap_fiber = interp_alphap(freq_array)
+
+            freq_diff = abs(freq_array - np.reshape(freq_array, (len(freq_array), 1)))
+            interp_cr = interp1d(raman_efficiency['frequency_offset'], raman_efficiency['cr'])
+            cr = interp_cr(freq_diff)
+
+            # z propagation axis
+            z = np.arange(0, fiber_length+1, z_resolution)
+
+            ode_function = lambda z, p: self._ode_stimulated_raman(z, p, alphap_fiber, freq_array, cr, prop_direct)
+            boundary_residual = lambda ya, yb: self._residuals_stimulated_raman(ya, yb, power_spectrum, prop_direct)
+            initial_guess_conditions = self._initial_guess_stimulated_raman(z, power_spectrum, alphap_fiber, prop_direct)
+
+            # ODE SOLVER
+            stimulated_raman_scattering = solve_bvp(ode_function, boundary_residual, z, initial_guess_conditions, tol=tolerance)
+
+            rho = (stimulated_raman_scattering.y.transpose() / power_spectrum).transpose()
+            rho = np.sqrt(rho)    # From power attenuation to field attenuation
+            setattr(stimulated_raman_scattering, 'frequency', freq_array)
+            setattr(stimulated_raman_scattering, 'z', stimulated_raman_scattering.x)
+            setattr(stimulated_raman_scattering, 'rho', rho)
+            setattr(stimulated_raman_scattering, 'power', stimulated_raman_scattering.y)
+            delattr(stimulated_raman_scattering, 'x')
+            delattr(stimulated_raman_scattering, 'y')
+
+            self.carriers = carriers
+            self.raman_pumps = raman_pumps
+            self._stimulated_raman_scattering = stimulated_raman_scattering
+
+        return self._stimulated_raman_scattering
+
+    def _residuals_stimulated_raman(self, ya, yb, power_spectrum, prop_direct):
+
+        computed_boundary_value = np.zeros(ya.size)
+
+        for index, direction in enumerate(prop_direct):
+            if direction == +1:
+                computed_boundary_value[index] = ya[index]
+            else:
+                computed_boundary_value[index] = yb[index]
+
+        return power_spectrum - computed_boundary_value
+
+    def _initial_guess_stimulated_raman(self, z, power_spectrum, alphap_fiber, prop_direct):
+        """ Computes the initial guess knowing the boundary conditions
+        :param z: patial axis [m]. numpy array
+        :param power_spectrum: power in each frequency slice [W].    Frequency axis is defined by freq_array. numpy array
+        :param alphap_fiber: frequency dependent fiber attenuation of signal power [1/m]. Frequency defined by freq_array. numpy array
+        :param prop_direct: indicates the propagation direction of each power slice in power_spectrum:
+        +1 for forward propagation and -1 for backward propagation. Frequency defined by freq_array. numpy array
+        :return: power_guess: guess on the initial conditions [W]. The first ndarray index identifies the frequency slice,
+        the second ndarray index identifies the step in z. ndarray
+        """
+
+        power_guess = np.empty((power_spectrum.size, z.size))
+        for f_index, power_slice in enumerate(power_spectrum):
+            if prop_direct[f_index] == +1:
+                power_guess[f_index, :] = np.exp(-alphap_fiber[f_index] * z) * power_slice
+            else:
+                power_guess[f_index, :] = np.exp(-alphap_fiber[f_index] * z[::-1]) * power_slice
+
+        return power_guess
+
+    def _ode_stimulated_raman(self, z, power_spectrum, alphap_fiber, freq_array, cr_raman_matrix, prop_direct):
+        """ Aim of ode_raman is to implement the set of ordinary differential equations (ODEs) describing the Raman effect.
+        :param z: spatial axis (unused).
+        :param power_spectrum: power in each frequency slice [W].    Frequency axis is defined by freq_array. numpy array. Size n
+        :param alphap_fiber: frequency dependent fiber attenuation of signal power [1/m]. Frequency defined by freq_array. numpy array. Size n
+        :param freq_array: reference frequency axis [Hz]. numpy array. Size n
+        :param cr_raman: Cr(f) Raman gain efficiency variation in frequency [1/W/m]. Frequency defined by freq_array. numpy ndarray. Size nxn
+        :param prop_direct: indicates the propagation direction of each power slice in power_spectrum:
+        +1 for forward propagation and -1 for backward propagation. Frequency defined by freq_array. numpy array. Size n
+        :return: dP/dz: the power variation in dz [W/m]. numpy array. Size n
+        """
+
+        dpdz = np.nan * np.ones(power_spectrum.shape)
+        for f_ind, power in enumerate(power_spectrum):
+            cr_raman = cr_raman_matrix[f_ind, :]
+            vibrational_loss = freq_array[f_ind] / freq_array[:f_ind]
+
+            for z_ind, power_sample in enumerate(power):
+                raman_gain = np.sum(cr_raman[f_ind+1:] * power_spectrum[f_ind+1:, z_ind])
+                raman_loss = np.sum(vibrational_loss * cr_raman[:f_ind] * power_spectrum[:f_ind, z_ind])
+
+                dpdz_element = prop_direct[f_ind] * (-alphap_fiber[f_ind] + raman_gain - raman_loss) * power_sample
+                dpdz[f_ind][z_ind] = dpdz_element
+
+        return np.vstack(dpdz)
+
+class NliSolver:
+    """ This class implements the NLI models.
+        Model and method can be specified in `self.nli_params.method`.
+        List of implemented methods:
+        'gn_model_analytic': brute force triple integral solution
+        'ggn_spectrally_separated_xpm_spm': XPM plus SPM
+    """
+
+    def __init__(self, nli_params=None, fiber_params=None):
+        """ Initialize the fiber object with its physical parameters
+        """
+        self.fiber_params = fiber_params
+        self.nli_params = nli_params
+        self.stimulated_raman_scattering = None
+
+    @property
+    def fiber_params(self):
+        return self._fiber_params
+
+    @fiber_params.setter
+    def fiber_params(self, fiber_params):
+        self._fiber_params = fiber_params
+
+    @property
+    def stimulated_raman_scattering(self):
+        return self._stimulated_raman_scattering
+
+    @stimulated_raman_scattering.setter
+    def stimulated_raman_scattering(self, stimulated_raman_scattering):
+        self._stimulated_raman_scattering = stimulated_raman_scattering
+
+    @property
+    def nli_params(self):
+        return self._nli_params
+
+    @nli_params.setter
+    def nli_params(self, nli_params):
+        """
+        :param model_params: namedtuple containing the parameters used to compute the NLI.
+        """
+        self._nli_params = nli_params
+
+    def compute_nli(self, carrier, *carriers):
+        """ Compute NLI power generated by the WDM comb `*carriers` on the channel under test `carrier`
+        at the end of the fiber span.
+        """
+        if 'gn_model_analytic' == self.nli_params.nli_method_name.lower():
+            carrier_nli = self._gn_analytic(carrier, *carriers)
+        elif 'ggn_spectrally_separated' in self.nli_params.nli_method_name.lower():
+            eta_matrix = self._compute_eta_matrix(carrier, *carriers)
+            carrier_nli = self._carrier_nli_from_eta_matrix(eta_matrix, carrier, *carriers)
+        else:
+            raise ValueError(f'Method {self.nli_params.method_nli} not implemented.')
+
+        return carrier_nli
+
+    @staticmethod
+    def _carrier_nli_from_eta_matrix(eta_matrix, carrier, *carriers):
+        carrier_nli = 0
+        for pump_carrier_1 in carriers:
+            for pump_carrier_2 in carriers:
+                carrier_nli += eta_matrix[pump_carrier_1.channel_number-1, pump_carrier_2.channel_number-1] * \
+                               pump_carrier_1.power.signal * pump_carrier_2.power.signal
+        carrier_nli *= carrier.power.signal
+
+        return carrier_nli
+
+    def _compute_eta_matrix(self, carrier_cut, *carriers):
+        cut_index = carrier_cut.channel_number - 1
+        # Matrix initialization
+        matrix_size = max(carriers, key=lambda x: getattr(x, 'channel_number')).channel_number
+        eta_matrix = np.zeros(shape=(matrix_size, matrix_size))
+
+        # SPM
+        logger.debug(f'Start computing SPM on channel #{carrier_cut.channel_number}')
+        # SPM GGN
+        if 'ggn' in self.nli_params.nli_method_name.lower():
+            partial_nli = self._generalized_spectrally_separated_spm(carrier_cut)
+        # SPM GN
+        elif 'gn' in self.nli_params.nli_method_name.lower():
+            partial_nli = self._gn_analytic(carrier_cut, *[carrier_cut])
+        eta_matrix[cut_index, cut_index] = partial_nli / (carrier_cut.power.signal**3)
+
+        # XPM
+        for pump_carrier in carriers:
+            pump_index = pump_carrier.channel_number - 1
+            if not (cut_index == pump_index):
+                logger.debug(f'Start computing XPM on channel #{carrier_cut.channel_number} '
+                             f'from channel #{pump_carrier.channel_number}')
+                # XPM GGN
+                if 'ggn' in self.nli_params.nli_method_name.lower():
+                    partial_nli = self._generalized_spectrally_separated_xpm(carrier_cut, pump_carrier)
+                # XPM GGN
+                elif 'gn' in self.nli_params.nli_method_name.lower():
+                    partial_nli = self._gn_analytic(carrier_cut, *[pump_carrier])
+                eta_matrix[pump_index, pump_index] = partial_nli /\
+                                                     (carrier_cut.power.signal * pump_carrier.power.signal**2)
+        return eta_matrix
+
+    # Methods for computing GN-model
+    def _gn_analytic(self, carrier, *carriers):
+        """ Computes the nonlinear interference power on a single carrier.
+        The method uses eq. 120 from arXiv:1209.0394.
+        :param carrier: the signal under analysis
+        :param carriers: the full WDM comb
+        :return: carrier_nli: the amount of nonlinear interference in W on the carrier under analysis
+        """
+        alpha = self.fiber_params.alpha0() / 2
+        beta2 = self.fiber_params.beta2
+        gamma = self.fiber_params.gamma
+        length = self.fiber_params.length
+        effective_length = (1 - np.exp(-2 * alpha * length)) / (2 * alpha)
+        asymptotic_length = 1 / (2 * alpha)
+
+        g_nli = 0
+        for interfering_carrier in carriers:
+            g_interfearing = interfering_carrier.power.signal / interfering_carrier.baud_rate
+            g_signal = carrier.power.signal / carrier.baud_rate
+            g_nli += g_interfearing**2 * g_signal \
+                * _psi(carrier, interfering_carrier, beta2=self.fiber_params.beta2, asymptotic_length=1/self.fiber_params.alpha0())
+        g_nli *= (16.0 / 27.0) * (gamma * effective_length)**2 /\
+                 (2 * np.pi * abs(beta2) * asymptotic_length)
+        carrier_nli = carrier.baud_rate * g_nli
+        return carrier_nli
+
+    # Methods for computing the GGN-model
+    def _generalized_spectrally_separated_spm(self, carrier):
+        f_cut_resolution = self.nli_params.f_cut_resolution['delta_0']
+        f_eval = carrier.frequency
+        g_cut = (carrier.power.signal / carrier.baud_rate)
+
+        spm_nli = carrier.baud_rate * (16.0 / 27.0) * self.fiber_params.gamma**2 * g_cut**3 * \
+                  self._generalized_psi(carrier, carrier, f_eval, f_cut_resolution, f_cut_resolution)
+        return spm_nli
+
+    def _generalized_spectrally_separated_xpm(self, carrier_cut, pump_carrier):
+        delta_index = pump_carrier.channel_number - carrier_cut.channel_number
+        f_cut_resolution = self.nli_params.f_cut_resolution[f'delta_{delta_index}']
+        f_pump_resolution = self.nli_params.f_pump_resolution
+        f_eval = carrier_cut.frequency
+        g_pump = (pump_carrier.power.signal / pump_carrier.baud_rate)
+        g_cut = (carrier_cut.power.signal / carrier_cut.baud_rate)
+        frequency_offset_threshold = self._frequency_offset_threshold(pump_carrier.baud_rate)
+        if abs(carrier_cut.frequency - pump_carrier.frequency) <= frequency_offset_threshold:
+            xpm_nli = carrier_cut.baud_rate * (16.0 / 27.0) * self.fiber_params.gamma**2 * g_pump**2 * g_cut * \
+                      2 * self._generalized_psi(carrier_cut, pump_carrier, f_eval, f_cut_resolution, f_pump_resolution)
+        else:
+            xpm_nli = carrier_cut.baud_rate * (16.0 / 27.0) * self.fiber_params.gamma**2 * g_pump**2 * g_cut * \
+                      2 * self._fast_generalized_psi(carrier_cut, pump_carrier, f_eval, f_cut_resolution)
+        return xpm_nli
+
+    def _fast_generalized_psi(self, carrier_cut, pump_carrier, f_eval, f_cut_resolution):
+        """ It computes the generalized psi function similarly to the one used in the GN model
+        :return: generalized_psi
+        """
+        # Fiber parameters
+        alpha0 = self.fiber_params.alpha0(f_eval)
+        beta2 = self.fiber_params.beta2
+        beta3 = self.fiber_params.beta3
+        f_ref_beta = self.fiber_params.f_ref_beta
+        z = self.stimulated_raman_scattering.z
+        frequency_rho = self.stimulated_raman_scattering.frequency
+        rho_norm = self.stimulated_raman_scattering.rho * np.exp(np.abs(alpha0) * z / 2)
+        if len(frequency_rho) == 1:
+            rho_function = lambda f: rho_norm[0, :]
+        else:
+            rho_function = interp1d(frequency_rho, rho_norm, axis=0, fill_value='extrapolate')
+        rho_norm_pump = rho_function(pump_carrier.frequency)
+
+        f1_array = np.array([pump_carrier.frequency - (pump_carrier.baud_rate * (1 + pump_carrier.roll_off) / 2),
+                             pump_carrier.frequency + (pump_carrier.baud_rate * (1 + pump_carrier.roll_off) / 2)])
+        f2_array = np.arange(carrier_cut.frequency,
+                             carrier_cut.frequency + (carrier_cut.baud_rate * (1 + carrier_cut.roll_off) / 2),
+                             f_cut_resolution)  # Only positive f2 is used since integrand_f2 is symmetric
+
+        integrand_f1 = np.zeros(len(f1_array))
+        for f1_index, f1 in enumerate(f1_array):
+            delta_beta = 4 * np.pi**2 * (f1 - f_eval) * (f2_array - f_eval) * \
+                         (beta2 + np.pi * beta3 * (f1 + f2_array - 2 * f_ref_beta))
+            integrand_f2 = self._generalized_rho_nli(delta_beta, rho_norm_pump, z, alpha0)
+            integrand_f1[f1_index] = 2 * np.trapz(integrand_f2, f2_array)  # 2x since integrand_f2 is symmetric in f2
+        generalized_psi = 0.5 * sum(integrand_f1) * pump_carrier.baud_rate
+        return generalized_psi
+
+    def _generalized_psi(self, carrier_cut, pump_carrier, f_eval, f_cut_resolution, f_pump_resolution):
+        """ It computes the generalized psi function similarly to the one used in the GN model
+        :return: generalized_psi
+        """
+        # Fiber parameters
+        alpha0 = self.fiber_params.alpha0(f_eval)
+        beta2 = self.fiber_params.beta2
+        beta3 = self.fiber_params.beta3
+        f_ref_beta = self.fiber_params.f_ref_beta
+        z = self.stimulated_raman_scattering.z
+        frequency_rho = self.stimulated_raman_scattering.frequency
+        rho_norm = self.stimulated_raman_scattering.rho * np.exp(np.abs(alpha0) * z / 2)
+        if len(frequency_rho) == 1:
+            rho_function = lambda f: rho_norm[0, :]
+        else:
+            rho_function = interp1d(frequency_rho, rho_norm, axis=0, fill_value='extrapolate')
+        rho_norm_pump = rho_function(pump_carrier.frequency)
+
+        f1_array = np.arange(pump_carrier.frequency - (pump_carrier.baud_rate * (1 + pump_carrier.roll_off) / 2),
+                             pump_carrier.frequency + (pump_carrier.baud_rate * (1 + pump_carrier.roll_off) / 2),
+                             f_pump_resolution)
+        f2_array = np.arange(carrier_cut.frequency - (carrier_cut.baud_rate * (1 + carrier_cut.roll_off) / 2),
+                             carrier_cut.frequency + (carrier_cut.baud_rate * (1 + carrier_cut.roll_off) / 2),
+                             f_cut_resolution)
+        psd1 = raised_cosine_comb(f1_array, pump_carrier) * (pump_carrier.baud_rate / pump_carrier.power.signal)
+
+        integrand_f1 = np.zeros(len(f1_array))
+        for f1_index, (f1, psd1_sample) in enumerate(zip(f1_array, psd1)):
+            f3_array = f1 + f2_array - f_eval
+            psd2 = raised_cosine_comb(f2_array, carrier_cut) * (carrier_cut.baud_rate / carrier_cut.power.signal)
+            psd3 = raised_cosine_comb(f3_array, pump_carrier) * (pump_carrier.baud_rate / pump_carrier.power.signal)
+            ggg = psd1_sample * psd2 * psd3
+
+            delta_beta = 4 * np.pi**2 * (f1 - f_eval) * (f2_array - f_eval) * \
+                         (beta2 + np.pi * beta3 * (f1 + f2_array - 2 * f_ref_beta))
+
+            integrand_f2 = ggg * self._generalized_rho_nli(delta_beta, rho_norm_pump, z, alpha0)
+            integrand_f1[f1_index] = np.trapz(integrand_f2, f2_array)
+        generalized_psi = np.trapz(integrand_f1, f1_array)
+        return generalized_psi
+
+    @staticmethod
+    def _generalized_rho_nli(delta_beta, rho_norm_pump, z, alpha0):
+        w = 1j * delta_beta - alpha0
+        generalized_rho_nli = (rho_norm_pump[-1]**2 * np.exp(w * z[-1]) - rho_norm_pump[0]**2 * np.exp(w * z[0])) / w
+        for z_ind in range(0, len(z) - 1):
+            derivative_rho = (rho_norm_pump[z_ind + 1]**2 - rho_norm_pump[z_ind]**2) / (z[z_ind + 1] - z[z_ind])
+            generalized_rho_nli -= derivative_rho * (np.exp(w * z[z_ind + 1]) - np.exp(w * z[z_ind])) / (w**2)
+        generalized_rho_nli = np.abs(generalized_rho_nli)**2
+        return generalized_rho_nli
+
+    def _frequency_offset_threshold(self, symbol_rate):
+        k_ref = 5
+        beta2_ref = 21.3e-27
+        delta_f_ref = 50e9
+        rs_ref = 32e9
+        freq_offset_th = ((k_ref * delta_f_ref) * rs_ref * beta2_ref) / (self.fiber_params.beta2 * symbol_rate)
+        return freq_offset_th
+
+def _psi(carrier, interfering_carrier, beta2, asymptotic_length):
+    """Calculates eq. 123 from `arXiv:1209.0394 <https://arxiv.org/abs/1209.0394>`__"""
+
+    if carrier.channel_number == interfering_carrier.channel_number: # SCI, SPM
+        psi = np.arcsinh(0.5 * np.pi**2 * asymptotic_length * abs(beta2) * carrier.baud_rate**2)
+    else: # XCI, XPM
+        delta_f = carrier.frequency - interfering_carrier.frequency
+        psi = np.arcsinh(np.pi**2 * asymptotic_length * abs(beta2) *
+                         carrier.baud_rate * (delta_f + 0.5 * interfering_carrier.baud_rate))
+        psi -= np.arcsinh(np.pi**2 * asymptotic_length * abs(beta2) *
+                          carrier.baud_rate * (delta_f - 0.5 * interfering_carrier.baud_rate))
+    return psi

gnpy/core/utils.py

@@ -11,8 +11,8 @@ This module contains utility functions that are used with gnpy.
 
 import json
 
-import numpy as np
 from csv import writer
+import numpy as np
 from numpy import pi, cos, sqrt, log10
 from scipy import constants
 
@@ -199,3 +199,43 @@ def rrc(ffs, baud_rate, alpha):
     p_inds = np.where(np.logical_and(np.abs(ffs) > 0, np.abs(ffs) < l_lim))
     hf[p_inds] = 1
     return sqrt(hf)
+
+def merge_amplifier_restrictions(dict1, dict2):
+    """Updates contents of dicts recursively
+
+    >>> d1 = {'params': {'restrictions': {'preamp_variety_list': [], 'booster_variety_list': []}}}
+    >>> d2 = {'params': {'target_pch_out_db': -20}}
+    >>> merge_amplifier_restrictions(d1, d2)
+    {'params': {'restrictions': {'preamp_variety_list': [], 'booster_variety_list': []}, 'target_pch_out_db': -20}}
+
+    >>> d3 = {'params': {'restrictions': {'preamp_variety_list': ['foo'], 'booster_variety_list': ['bar']}}}
+    >>> merge_amplifier_restrictions(d1, d3)
+    {'params': {'restrictions': {'preamp_variety_list': [], 'booster_variety_list': []}}}
+    """
+
+    copy_dict1 = dict1.copy()
+    for key in dict2:
+        if key in dict1:
+            if isinstance(dict1[key], dict):
+                copy_dict1[key] = merge_amplifier_restrictions(copy_dict1[key], dict2[key])
+        else:
+            copy_dict1[key] = dict2[key]
+    return copy_dict1
+
+def silent_remove(this_list, elem):
+    """Remove matching elements from a list without raising ValueError
+
+    >>> li = [0, 1]
+    >>> li = silent_remove(li, 1)
+    >>> li
+    [0]
+    >>> li = silent_remove(li, 1)
+    >>> li
+    [0]
+    """
+
+    try:
+        this_list.remove(elem)
+    except ValueError:
+        pass
+    return this_list

json_structure_description.rst

@@ -165,6 +165,52 @@ Fiber element with its parameters:
             }
         ]
 
+RamanFiber element
+******************
+
+A special variant of the regular ``Fiber`` where the simulation engine accounts for the Raman effect.
+The newly added parameters are nested in the ``raman_efficiency`` dictionary.
+Its shape corresponds to typical properties of silica.
+More details are available from :cite:`curri_merit_2016`.
+
+The ``cr`` property is the normailzed Raman efficiency, so it is is (almost) independent of the fiber type, while the coefficient actually giving Raman gain is g_R=C_R/Aeff.
+
+The ``frequency_offset`` represents the spectral difference between the pumping photon and the one receiving energy.
+
+.. code-block:: json-object
+
+    "RamanFiber":[{
+      "type_variety": "SSMF",
+      "dispersion": 1.67e-05,
+      "gamma": 0.00127,
+      "raman_efficiency": {
+        "cr":[
+            0, 9.4E-06, 2.92E-05, 4.88E-05, 6.82E-05, 8.31E-05, 9.4E-05, 0.0001014, 0.0001069, 0.0001119,
+            0.0001217, 0.0001268, 0.0001365, 0.000149, 0.000165, 0.000181, 0.0001977, 0.0002192, 0.0002469,
+            0.0002749, 0.0002999, 0.0003206, 0.0003405, 0.0003592, 0.000374, 0.0003826, 0.0003841, 0.0003826,
+            0.0003802, 0.0003756, 0.0003549, 0.0003795, 0.000344, 0.0002933, 0.0002024, 0.0001158, 8.46E-05,
+            7.14E-05, 6.86E-05, 8.5E-05, 8.93E-05, 9.01E-05, 8.15E-05, 6.67E-05, 4.37E-05, 3.28E-05, 2.96E-05,
+            2.65E-05, 2.57E-05, 2.81E-05, 3.08E-05, 3.67E-05, 5.85E-05, 6.63E-05, 6.36E-05, 5.5E-05, 4.06E-05,
+            2.77E-05, 2.42E-05, 1.87E-05, 1.6E-05, 1.4E-05, 1.13E-05, 1.05E-05, 9.8E-06, 9.8E-06, 1.13E-05,
+            1.64E-05, 1.95E-05, 2.38E-05, 2.26E-05, 2.03E-05, 1.48E-05, 1.09E-05, 9.8E-06, 1.05E-05, 1.17E-05,
+            1.25E-05, 1.21E-05, 1.09E-05, 9.8E-06, 8.2E-06, 6.6E-06, 4.7E-06, 2.7E-06, 1.9E-06, 1.2E-06, 4E-07,
+            2E-07, 1E-07
+        ],
+        "frequency_offset":[
+          0, 0.5e12, 1e12, 1.5e12, 2e12, 2.5e12, 3e12, 3.5e12, 4e12, 4.5e12, 5e12, 5.5e12, 6e12, 6.5e12, 7e12,
+          7.5e12, 8e12, 8.5e12, 9e12, 9.5e12, 10e12, 10.5e12, 11e12, 11.5e12, 12e12, 12.5e12, 12.75e12,
+          13e12, 13.25e12, 13.5e12, 14e12, 14.5e12, 14.75e12, 15e12, 15.5e12, 16e12, 16.5e12, 17e12,
+          17.5e12, 18e12, 18.25e12, 18.5e12, 18.75e12, 19e12, 19.5e12, 20e12, 20.5e12, 21e12, 21.5e12,
+          22e12, 22.5e12, 23e12, 23.5e12, 24e12, 24.5e12, 25e12, 25.5e12, 26e12, 26.5e12, 27e12, 27.5e12, 28e12,
+          28.5e12, 29e12, 29.5e12, 30e12, 30.5e12, 31e12, 31.5e12, 32e12, 32.5e12, 33e12, 33.5e12, 34e12, 34.5e12,
+          35e12, 35.5e12, 36e12, 36.5e12, 37e12, 37.5e12, 38e12, 38.5e12, 39e12, 39.5e12, 40e12, 40.5e12, 41e12,
+          41.5e12, 42e12
+        ]
+        }
+      }
+    ]
+
+
 1.2.3 Roadm element
 *******************
 
@@ -227,6 +273,8 @@ Spectral information with its parameters:
 Transceiver element with its parameters. **”mode”** can contain multiple
 Transceiver operation formats.
 
+Note that ``OSNR`` parameter refers to the receiver's minimal OSNR threshold for a given mode.
+
 .. code-block:: json-object
 
     "Transceiver":[{
@@ -406,8 +454,51 @@ Fiber element with its parameters.
                 }
     }
 
+2.2.5. RamanFiber element
+*************************
 
-2.2.5. EDFA element
+.. code-block:: json
+
+    {
+      "uid": "Span1",
+      "type": "RamanFiber",
+      "type_variety": "SSMF",
+      "operational": {
+        "temperature": 283,
+        "raman_pumps": [
+          {
+            "power": 200e-3,
+            "frequency": 205e12,
+            "propagation_direction": "counterprop"
+          },
+          {
+            "power": 206e-3,
+            "frequency": 201e12,
+            "propagation_direction": "counterprop"
+          }
+        ]
+      },
+      "params": {
+        "type_variety": "SSMF",
+        "length": 80.0,
+        "loss_coef": 0.2,
+        "length_units": "km",
+        "att_in": 0,
+        "con_in": 0.5,
+        "con_out": 0.5
+      },
+      "metadata": {
+        "location": {
+          "latitude": 1,
+          "longitude": 0,
+          "city": null,
+          "region": ""
+        }
+      }
+    }
+
+
+2.2.6. EDFA element
 ********************
 
 EDFA element with its parameters.
@@ -443,3 +534,26 @@ corresponding to element **”uid”**
     {"from_node": "roadm Site_C",
     "to_node": "trx Site_C"
     }
+
+************************
+3. Simulation Parameters
+************************
+
+Additional details of the simulation are controlled via ``sim_params.json``:
+
+.. code-block:: json
+
+  {
+    "raman_computed_channels": [1, 18, 37, 56, 75],
+    "raman_parameters": {
+      "flag_raman": true,
+      "space_resolution": 10e3,
+      "tolerance": 1e-8
+    },
+    "nli_parameters": {
+      "nli_method_name": "ggn_spectrally_separated",
+      "wdm_grid_size": 50e9,
+      "dispersion_tolerance": 1,
+      "phase_shift_tollerance": 0.1
+    }
+  }

pytest.ini

@@ -1,2 +1,2 @@
 [pytest]
-addopts = -p no:warnings
+addopts = --doctest-modules

requirements.txt

@@ -1,44 +1,10 @@
-alabaster==0.7.12
+alabaster>=0.7.12,<1
-appdirs==1.4.3
+matplotlib>=3.1.0,<4
-atomicwrites==1.2.1
+networkx>=2.3,<3
-attrs==18.2.0
+numpy>=1.16.1,<2
-Babel==2.6.0
+Pygments>=2.4.2,<3
-black==18.9b0
+pytest>=4.0.0,<5
-certifi==2018.10.15
+scipy>=1.3.0,<2
-chardet==3.0.4
+Sphinx>=2.1.1,<3
-Click==7.0
+sphinxcontrib-bibtex>=0.4.2,<1
-cycler==0.10.0
+xlrd>=1.2.0,<2
-decorator==4.3.0
-docutils==0.14
-idna==2.7
-imagesize==1.1.0
-Jinja2==2.10
-kiwisolver==1.0.1
-latexcodec==1.0.5
-MarkupSafe==1.0
-matplotlib==3.0.0
-more-itertools==4.3.0
-networkx==2.2
-numpy==1.15.2
-oset==0.1.3
-packaging==18.0
-pluggy==0.7.1
-py==1.7.0
-pybtex==0.21
-pybtex-docutils==0.2.1
-Pygments==2.2.0
-pyparsing==2.2.2
-pytest==3.8.2
-python-dateutil==2.7.3
-pytz==2018.5
-PyYAML==3.13
-requests==2.19.1
-scipy==1.1.0
-six==1.11.0
-snowballstemmer==1.2.1
-Sphinx==1.8.1
-sphinxcontrib-bibtex==0.4.0
-sphinxcontrib-websupport==1.1.0
-toml==0.10.0
-urllib3==1.23
-xlrd==1.1.0

tests/LinkforTest.json

@@ -205,6 +205,36 @@
                     "longitude": 0
                 }
             }
+        },
+        {
+        "uid": "Att_B",
+        "type": "Fused",
+        "params":{
+            "loss":16
+        },
+        "metadata": {
+            "location": {
+              "latitude": 2.0,
+              "longitude": 1.0,
+              "city": "Corlay",
+              "region": "RLD"
+                }
+            }
+        },
+        {
+        "uid": "Att_F",
+        "type": "Fused",
+        "params":{
+            "loss":16
+        },
+        "metadata": {
+            "location": {
+              "latitude": 2.0,
+              "longitude": 1.0,
+              "city": "Corlay",
+              "region": "RLD"
+                }
+            }
         }
 
   ],
@@ -247,6 +277,10 @@
     },
     {
       "from_node": "Edfa5",
+      "to_node": "Att_F"
+    },
+    {
+      "from_node": "Att_F",
       "to_node": "trx F"
     },
     {
@@ -255,6 +289,10 @@
     },
     {
       "from_node": "Edfa1",
+      "to_node": "Att_B"
+    },
+    {
+      "from_node": "Att_B",
       "to_node": "trx B"
     }
   ]

tests/data/eqpt_config.json

@@ -49,7 +49,16 @@
             "p_max": 21,
             "nf0": 5,
             "allowed_for_design": true
-            }
+            },
+            {
+            "type_variety": "std_booster",
+            "type_def": "fixed_gain",
+            "gain_flatmax": 21,
+            "gain_min": 20,
+            "p_max": 21,
+            "nf0": 5,
+            "allowed_for_design": false
+            }            
       ],
       "Fiber":[{
             "type_variety": "SSMF",
@@ -75,8 +84,8 @@
             "target_pch_out_db": -20,
             "add_drop_osnr": 38,
             "restrictions": {
-                            "preamp_variety_list":["low_gain_preamp", "high_gain_preamp"],
+                            "preamp_variety_list":[],
-                            "booster_variety_list":["std_booster"]
+                            "booster_variety_list":[]
                             }    
             }],
       "SI":[{

tests/data/testTopology_expected.json

@@ -42,7 +42,7 @@
         "location": {
           "city": "Rennes_STA",
           "region": "RLD",
-          "latitude": 0.0,
+          "latitude": 4.0,
           "longitude": 0.0
         }
       },
@@ -54,7 +54,7 @@
         "location": {
           "city": "Brest_KLA",
           "region": "RLD",
-          "latitude": 4.0,
+          "latitude": 0.0,
           "longitude": 0.0
         }
       },
@@ -66,8 +66,8 @@
         "location": {
           "city": "a",
           "region": "",
-          "latitude": 0,
+          "latitude": 6.0,
-          "longitude": 0
+          "longitude": 0.0
         }
       },
       "type": "Transceiver"
@@ -78,8 +78,8 @@
         "location": {
           "city": "b",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 0.0
         }
       },
       "type": "Transceiver"
@@ -90,8 +90,8 @@
         "location": {
           "city": "c",
           "region": "",
-          "latitude": 0,
+          "latitude": 6.0,
-          "longitude": 0
+          "longitude": 1.0
         }
       },
       "type": "Transceiver"
@@ -102,8 +102,8 @@
         "location": {
           "city": "d",
           "region": "",
-          "latitude": 0,
+          "latitude": 6.0,
-          "longitude": 0
+          "longitude": 4.0
         }
       },
       "type": "Transceiver"
@@ -114,8 +114,8 @@
         "location": {
           "city": "e",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 4.0
         }
       },
       "type": "Transceiver"
@@ -126,8 +126,8 @@
         "location": {
           "city": "f",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 1.0
         }
       },
       "type": "Transceiver"
@@ -138,8 +138,8 @@
         "location": {
           "city": "g",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 3.0
         }
       },
       "type": "Transceiver"
@@ -150,8 +150,8 @@
         "location": {
           "city": "h",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 2.0
         }
       },
       "type": "Transceiver"
@@ -198,7 +198,7 @@
         "location": {
           "city": "Rennes_STA",
           "region": "RLD",
-          "latitude": 0.0,
+          "latitude": 4.0,
           "longitude": 0.0
         }
       },
@@ -210,7 +210,7 @@
         "location": {
           "city": "Brest_KLA",
           "region": "RLD",
-          "latitude": 4.0,
+          "latitude": 0.0,
           "longitude": 0.0
         }
       },
@@ -218,24 +218,40 @@
     },
     {
       "uid": "roadm a",
+      "params": {
+        "restrictions": {
+          "preamp_variety_list": [],
+          "booster_variety_list": [
+            "std_booster"
+          ]
+        }
+      },
       "metadata": {
         "location": {
           "city": "a",
           "region": "",
-          "latitude": 0,
+          "latitude": 6.0,
-          "longitude": 0
+          "longitude": 0.0
         }
       },
       "type": "Roadm"
     },
     {
       "uid": "roadm b",
+      "params": {
+        "restrictions": {
+          "preamp_variety_list": [
+            "std_low_gain"
+          ],
+          "booster_variety_list": []
+        }
+      },
       "metadata": {
         "location": {
           "city": "b",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 0.0
         }
       },
       "type": "Roadm"
@@ -246,8 +262,8 @@
         "location": {
           "city": "c",
           "region": "",
-          "latitude": 0,
+          "latitude": 6.0,
-          "longitude": 0
+          "longitude": 1.0
         }
       },
       "type": "Roadm"
@@ -258,8 +274,8 @@
         "location": {
           "city": "d",
           "region": "",
-          "latitude": 0,
+          "latitude": 6.0,
-          "longitude": 0
+          "longitude": 4.0
         }
       },
       "type": "Roadm"
@@ -270,8 +286,8 @@
         "location": {
           "city": "e",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 4.0
         }
       },
       "type": "Roadm"
@@ -282,8 +298,8 @@
         "location": {
           "city": "f",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 1.0
         }
       },
       "type": "Roadm"
@@ -294,8 +310,8 @@
         "location": {
           "city": "g",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 3.0
         }
       },
       "type": "Roadm"
@@ -306,8 +322,8 @@
         "location": {
           "city": "h",
           "region": "",
-          "latitude": 0,
+          "latitude": 5.0,
-          "longitude": 0
+          "longitude": 2.0
         }
       },
       "type": "Roadm"
@@ -342,7 +358,7 @@
         "location": {
           "city": "Morlaix",
           "region": "RLD",
-          "latitude": 3.0,
+          "latitude": 1.0,
           "longitude": 0.0
         }
       },
@@ -378,7 +394,7 @@
         "location": {
           "city": "Morlaix",
           "region": "RLD",
-          "latitude": 3.0,
+          "latitude": 1.0,
           "longitude": 0.0
         }
       },
@@ -460,7 +476,7 @@
       "uid": "fiber (Lannion_CAS → Stbrieuc)-F056",
       "metadata": {
         "location": {
-          "latitude": 1.5,
+          "latitude": 2.5,
           "longitude": 0.0
         }
       },
@@ -478,7 +494,7 @@
       "uid": "fiber (Stbrieuc → Rennes_STA)-F057",
       "metadata": {
         "location": {
-          "latitude": 0.5,
+          "latitude": 3.5,
           "longitude": 0.0
         }
       },
@@ -496,7 +512,7 @@
       "uid": "fiber (Lannion_CAS → Morlaix)-F059",
       "metadata": {
         "location": {
-          "latitude": 2.5,
+          "latitude": 1.5,
           "longitude": 0.0
         }
       },
@@ -514,7 +530,7 @@
       "uid": "fiber (Morlaix → Brest_KLA)-F060",
       "metadata": {
         "location": {
-          "latitude": 3.5,
+          "latitude": 0.5,
           "longitude": 0.0
         }
       },
@@ -532,8 +548,8 @@
       "uid": "fiber (Brest_KLA → Quimper)-",
       "metadata": {
         "location": {
-          "latitude": 2.5,
+          "latitude": 0.0,
-          "longitude": 0.5
+          "longitude": 1.5
         }
       },
       "type": "Fiber",
@@ -550,8 +566,8 @@
       "uid": "fiber (Quimper → Lorient_KMA)-",
       "metadata": {
         "location": {
-          "latitude": 1.5,
+          "latitude": 1.0,
-          "longitude": 2.0
+          "longitude": 3.0
         }
       },
       "type": "Fiber",
@@ -568,8 +584,8 @@
       "uid": "fiber (Ploermel → Vannes_KBE)-",
       "metadata": {
         "location": {
-          "latitude": 1.5,
+          "latitude": 3.0,
-          "longitude": 3.0
+          "longitude": 4.0
         }
       },
       "type": "Fiber",
@@ -586,8 +602,8 @@
       "uid": "fiber (Ploermel → Rennes_STA)-",
       "metadata": {
         "location": {
-          "latitude": 0.5,
+          "latitude": 4.0,
-          "longitude": 1.0
+          "longitude": 2.0
         }
       },
       "type": "Fiber",
@@ -604,7 +620,7 @@
       "uid": "fiber (a → b)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.5,
           "longitude": 0.0
         }
       },
@@ -622,8 +638,8 @@
       "uid": "fiber (a → c)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 6.0,
-          "longitude": 0.0
+          "longitude": 0.5
         }
       },
       "type": "Fiber",
@@ -640,14 +656,14 @@
       "uid": "fiber (c → d)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 6.0,
-          "longitude": 0.0
+          "longitude": 2.5
         }
       },
       "type": "Fiber",
       "type_variety": "SSMF",
       "params": {
-        "length": 50.0,
+        "length": 10.0,
         "length_units": "km",
         "loss_coef": 0.2,
         "con_in": null,
@@ -658,8 +674,8 @@
       "uid": "fiber (c → f)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.5,
-          "longitude": 0.0
+          "longitude": 1.0
         }
       },
       "type": "Fiber",
@@ -676,8 +692,8 @@
       "uid": "fiber (b → f)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.0,
-          "longitude": 0.0
+          "longitude": 0.5
         }
       },
       "type": "Fiber",
@@ -694,8 +710,8 @@
       "uid": "fiber (e → d)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.5,
-          "longitude": 0.0
+          "longitude": 4.0
         }
       },
       "type": "Fiber",
@@ -712,8 +728,8 @@
       "uid": "fiber (e → g)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.0,
-          "longitude": 0.0
+          "longitude": 3.5
         }
       },
       "type": "Fiber",
@@ -730,8 +746,8 @@
       "uid": "fiber (f → h)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.0,
-          "longitude": 0.0
+          "longitude": 1.5
         }
       },
       "type": "Fiber",
@@ -748,8 +764,8 @@
       "uid": "fiber (h → g)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.0,
-          "longitude": 0.0
+          "longitude": 2.5
         }
       },
       "type": "Fiber",
@@ -838,7 +854,7 @@
       "uid": "fiber (Stbrieuc → Lannion_CAS)-F056",
       "metadata": {
         "location": {
-          "latitude": 1.5,
+          "latitude": 2.5,
           "longitude": 0.0
         }
       },
@@ -856,7 +872,7 @@
       "uid": "fiber (Rennes_STA → Stbrieuc)-F057",
       "metadata": {
         "location": {
-          "latitude": 0.5,
+          "latitude": 3.5,
           "longitude": 0.0
         }
       },
@@ -874,7 +890,7 @@
       "uid": "fiber (Morlaix → Lannion_CAS)-F059",
       "metadata": {
         "location": {
-          "latitude": 2.5,
+          "latitude": 1.5,
           "longitude": 0.0
         }
       },
@@ -892,7 +908,7 @@
       "uid": "fiber (Brest_KLA → Morlaix)-F060",
       "metadata": {
         "location": {
-          "latitude": 3.5,
+          "latitude": 0.5,
           "longitude": 0.0
         }
       },
@@ -910,8 +926,8 @@
       "uid": "fiber (Quimper → Brest_KLA)-",
       "metadata": {
         "location": {
-          "latitude": 2.5,
+          "latitude": 0.0,
-          "longitude": 0.5
+          "longitude": 1.5
         }
       },
       "type": "Fiber",
@@ -928,8 +944,8 @@
       "uid": "fiber (Lorient_KMA → Quimper)-",
       "metadata": {
         "location": {
-          "latitude": 1.5,
+          "latitude": 1.0,
-          "longitude": 2.0
+          "longitude": 3.0
         }
       },
       "type": "Fiber",
@@ -946,8 +962,8 @@
       "uid": "fiber (Vannes_KBE → Ploermel)-",
       "metadata": {
         "location": {
-          "latitude": 1.5,
+          "latitude": 3.0,
-          "longitude": 3.0
+          "longitude": 4.0
         }
       },
       "type": "Fiber",
@@ -964,8 +980,8 @@
       "uid": "fiber (Rennes_STA → Ploermel)-",
       "metadata": {
         "location": {
-          "latitude": 0.5,
+          "latitude": 4.0,
-          "longitude": 1.0
+          "longitude": 2.0
         }
       },
       "type": "Fiber",
@@ -982,7 +998,7 @@
       "uid": "fiber (b → a)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.5,
           "longitude": 0.0
         }
       },
@@ -1000,8 +1016,8 @@
       "uid": "fiber (c → a)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 6.0,
-          "longitude": 0.0
+          "longitude": 0.5
         }
       },
       "type": "Fiber",
@@ -1018,14 +1034,14 @@
       "uid": "fiber (d → c)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 6.0,
-          "longitude": 0.0
+          "longitude": 2.5
         }
       },
       "type": "Fiber",
       "type_variety": "SSMF",
       "params": {
-        "length": 50.0,
+        "length": 10.0,
         "length_units": "km",
         "loss_coef": 0.2,
         "con_in": null,
@@ -1036,8 +1052,8 @@
       "uid": "fiber (f → c)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.5,
-          "longitude": 0.0
+          "longitude": 1.0
         }
       },
       "type": "Fiber",
@@ -1054,8 +1070,8 @@
       "uid": "fiber (f → b)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.0,
-          "longitude": 0.0
+          "longitude": 0.5
         }
       },
       "type": "Fiber",
@@ -1072,8 +1088,8 @@
       "uid": "fiber (d → e)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.5,
-          "longitude": 0.0
+          "longitude": 4.0
         }
       },
       "type": "Fiber",
@@ -1090,8 +1106,8 @@
       "uid": "fiber (g → e)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.0,
-          "longitude": 0.0
+          "longitude": 3.5
         }
       },
       "type": "Fiber",
@@ -1108,8 +1124,8 @@
       "uid": "fiber (h → f)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.0,
-          "longitude": 0.0
+          "longitude": 1.5
         }
       },
       "type": "Fiber",
@@ -1126,8 +1142,8 @@
       "uid": "fiber (g → h)-",
       "metadata": {
         "location": {
-          "latitude": 0.0,
+          "latitude": 5.0,
-          "longitude": 0.0
+          "longitude": 2.5
         }
       },
       "type": "Fiber",
@@ -1178,25 +1194,6 @@
         "out_voa": null
       }
     },
-    {
-      "uid": "east edfa in Stbrieuc to Rennes_STA",
-      "metadata": {
-        "location": {
-          "city": "Stbrieuc",
-          "region": "RLD",
-          "latitude": 1.0,
-          "longitude": 0.0
-        }
-      },
-      "type": "Edfa",
-      "type_variety": "std_medium_gain",
-      "operational": {
-        "gain_target": 18.5,
-        "delta_p": null,
-        "tilt_target": 0,
-        "out_voa": null
-      }
-    },
     {
       "uid": "east edfa in Lannion_CAS to Morlaix",
       "metadata": {
@@ -1216,13 +1213,32 @@
         "out_voa": 0.5
       }
     },
+    {
+      "uid": "east edfa in Stbrieuc to Rennes_STA",
+      "metadata": {
+        "location": {
+          "city": "Stbrieuc",
+          "region": "RLD",
+          "latitude": 3.0,
+          "longitude": 0.0
+        }
+      },
+      "type": "Edfa",
+      "type_variety": "std_medium_gain",
+      "operational": {
+        "gain_target": 18.5,
+        "delta_p": null,
+        "tilt_target": 0,
+        "out_voa": null
+      }
+    },
     {
       "uid": "east edfa in Brest_KLA to Quimper",
       "metadata": {
         "location": {
           "city": "Brest_KLA",
           "region": "RLD",
-          "latitude": 4.0,
+          "latitude": 0.0,
           "longitude": 0.0
         }
       },
@@ -1241,8 +1257,8 @@
         "location": {
           "city": "Ploermel",
           "region": "RLD",
-          "latitude": 1.0,
+          "latitude": 4.0,
-          "longitude": 2.0
+          "longitude": 4.0
         }
       },
       "type": "Edfa",
@@ -1293,22 +1309,18 @@
       }
     },
     {
-      "uid": "west edfa in Quimper to Lorient_KMA",
+      "uid": "east edfa in c to d",
       "metadata": {
         "location": {
-          "city": "Quimper",
+          "city": "c",
-          "region": "RLD",
+          "region": "",
-          "latitude": 1.0,
+          "latitude": 6.0,
           "longitude": 1.0
         }
       },
-      "type": "Edfa",
+      "type": "Fused",
-      "type_variety": "std_low_gain",
+      "params": {
-      "operational": {
+        "loss": 0
-        "gain_target": 19.0,
-        "delta_p": null,
-        "tilt_target": 0,
-        "out_voa": null
       }
     }
   ],
@@ -1499,10 +1511,6 @@
     },
     {
       "from_node": "fiber (Lorient_KMA → Quimper)-",
-      "to_node": "west edfa in Quimper to Lorient_KMA"
-    },
-    {
-      "from_node": "west edfa in Quimper to Lorient_KMA",
       "to_node": "fiber (Quimper → Brest_KLA)-"
     },
     {
@@ -1559,6 +1567,10 @@
     },
     {
       "from_node": "roadm c",
+      "to_node": "east edfa in c to d"
+    },
+    {
+      "from_node": "east edfa in c to d",
       "to_node": "fiber (c → d)-"
     },
     {
@@ -1766,4 +1778,4 @@
       "to_node": "trx h"
     }
   ]
-}
+}

tests/data/test_network.json

@@ -77,7 +77,22 @@
                     "longitude": 0
                 }
             }
-        },   
+        },
+        {
+        "uid": "Att_B",
+        "type": "Fused",
+        "params":{
+            "loss":16
+        }, 
+        "metadata": {
+            "location": {
+              "latitude": 2.0,
+              "longitude": 1.0,
+              "city": "Corlay",
+              "region": "RLD"
+                }
+            }
+        },
         {
             "uid": "Site_B",
             "type": "Transceiver",
@@ -110,6 +125,10 @@
         },        
         {
             "from_node": "Edfa2",
+            "to_node": "Att_B"
+        },
+        {
+            "from_node": "Att_B",
             "to_node": "Site_B"
         }
 

tests/test_automaticmodefeature.py

@@ -72,7 +72,7 @@ def test_automaticmodefeature(net,eqpt,serv,expected_mode):
         if pathreq.baud_rate is not None:
             print(pathreq.format)
             path_res_list.append(pathreq.format)
-            total_path = propagate(total_path,pathreq,equipment, show=False)
+            total_path = propagate(total_path,pathreq,equipment)
         else:
             total_path,mode = propagate_and_optimize_mode(total_path,pathreq,equipment)
             # if no baudrate satisfies spacing, no mode is returned and an empty path is returned

tests/test_propagation.py

@@ -18,6 +18,8 @@ from numpy import mean
 
 #network_file_name = 'tests/test_network.json'
 network_file_name = Path(__file__).parent.parent / 'tests/LinkforTest.json'
+#TODO: note that this json entries has a weird topology since EDfa1 has a possible branch on a receiver B
+# this might not pass future tests/ code updates
 #network_file_name = Path(__file__).parent.parent / 'examples/edfa_example_network.json'
 eqpt_library_name = Path(__file__).parent.parent / 'tests/data/eqpt_config.json'
 

tests/test_roadm_restrictions.py

@@ -0,0 +1,203 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# @Author: Esther Le Rouzic
+# @Date:   2019-05-22
+"""
+@author: esther.lerouzic
+checks that fused placed in amp type is correctly converted to a fused element instead of an edfa
+and that no additional amp is added.
+checks that restrictions in roadms are correctly applied during autodesign
+
+"""
+
+from pathlib import Path
+import pytest
+from gnpy.core.utils import lin2db, load_json
+from gnpy.core.elements import Fused, Roadm, Edfa
+from gnpy.core.equipment import load_equipment, Amp, automatic_nch
+from gnpy.core.network import network_from_json, build_network
+
+
+TEST_DIR = Path(__file__).parent
+EQPT_LIBRARY_NAME = TEST_DIR / 'data/eqpt_config.json'
+NETWORK_FILE_NAME = TEST_DIR / 'data/testTopology_expected.json'
+# adding tests to check the roadm restrictions
+
+# mark node_uid amps as fused for testing purpose
+@pytest.mark.parametrize("node_uid", ['east edfa in Lannion_CAS to Stbrieuc'])
+def test_no_amp_feature(node_uid):
+    ''' Check that booster is not placed on a roadm if fused is specified
+        test_parser covers partly this behaviour. This test should guaranty that the
+        feature is preserved even if convert is changed
+    '''
+    equipment = load_equipment(EQPT_LIBRARY_NAME)
+    json_network = load_json(NETWORK_FILE_NAME)
+
+    for elem in json_network['elements']:
+        if elem['uid'] == node_uid:
+            #replace edfa node by a fused node in the topology
+            elem['type'] = 'Fused'
+            elem.pop('type_variety')
+            elem.pop('operational')
+            elem['params'] = {'loss': 0}
+
+            next_node_uid = next(conn['to_node'] for conn in json_network['connections'] \
+                                 if conn['from_node'] == node_uid)
+            previous_node_uid = next(conn['from_node'] for conn in json_network['connections'] \
+                                 if conn['to_node'] == node_uid)
+
+    network = network_from_json(json_network, equipment)
+    # Build the network once using the default power defined in SI in eqpt config
+    # power density : db2linp(ower_dbm": 0)/power_dbm": 0 * nb channels as defined by
+    # spacing, f_min and f_max
+    p_db = equipment['SI']['default'].power_dbm
+    p_total_db = p_db + lin2db(automatic_nch(equipment['SI']['default'].f_min,\
+        equipment['SI']['default'].f_max, equipment['SI']['default'].spacing))
+
+    build_network(network, equipment, p_db, p_total_db)
+
+    node = next(nd for nd in network.nodes() if nd.uid == node_uid)
+    next_node = next(network.successors(node))
+    previous_node = next(network.predecessors(node))
+
+    if not isinstance(node, Fused):
+        raise AssertionError()
+    if not node.params.loss == 0.0:
+        raise AssertionError()
+    if not next_node_uid == next_node.uid:
+        raise AssertionError()
+    if not previous_node_uid == previous_node.uid:
+        raise AssertionError()
+
+@pytest.fixture()
+def equipment():
+    """init transceiver class to access snr and osnr calculations"""
+    equipment = load_equipment(EQPT_LIBRARY_NAME)
+    # define some booster and preamps
+    restrictions_list = [
+        {
+            'type_variety': 'booster_medium_gain',
+            'type_def': 'variable_gain',
+            'gain_flatmax': 25,
+            'gain_min': 15,
+            'p_max': 21,
+            'nf_min': 5.8,
+            'nf_max': 10,
+            'out_voa_auto': False,
+            'allowed_for_design': False
+        },
+        {
+            'type_variety': 'preamp_medium_gain',
+            'type_def': 'variable_gain',
+            'gain_flatmax': 26,
+            'gain_min': 15,
+            'p_max': 23,
+            'nf_min': 6,
+            'nf_max': 10,
+            'out_voa_auto': False,
+            'allowed_for_design': False
+        },
+        {
+            'type_variety': 'preamp_high_gain',
+            'type_def': 'variable_gain',
+            'gain_flatmax': 35,
+            'gain_min': 25,
+            'p_max': 21,
+            'nf_min': 5.5,
+            'nf_max': 7,
+            'out_voa_auto': False,
+            'allowed_for_design': False
+        },
+        {
+            'type_variety': 'preamp_low_gain',
+            'type_def': 'variable_gain',
+            'gain_flatmax': 16,
+            'gain_min': 8,
+            'p_max': 23,
+            'nf_min': 6.5,
+            'nf_max': 11,
+            'out_voa_auto': False,
+            'allowed_for_design': False
+        }]
+    # add them to the library
+    for entry in restrictions_list:
+        equipment['Edfa'][entry['type_variety']] = Amp.from_json(EQPT_LIBRARY_NAME, **entry)
+    return equipment
+
+
+@pytest.mark.parametrize("restrictions", [
+    {
+        'preamp_variety_list':[],
+        'booster_variety_list':[]
+    },
+    {
+        'preamp_variety_list':[],
+        'booster_variety_list':['booster_medium_gain']
+    },
+    {
+        'preamp_variety_list':['preamp_medium_gain', 'preamp_high_gain', 'preamp_low_gain'],
+        'booster_variety_list':[]
+    }])
+def test_restrictions(restrictions, equipment):
+    ''' test that restriction is correctly applied if provided in eqpt_config and if no Edfa type
+    were provided in the network json
+    '''
+    # add restrictions
+    equipment['Roadm']['default'].restrictions = restrictions
+    # build network
+    json_network = load_json(NETWORK_FILE_NAME)
+    network = network_from_json(json_network, equipment)
+
+    amp_nodes_nobuild_uid = [nd.uid for nd in network.nodes() \
+        if isinstance(nd, Edfa) and isinstance(next(network.predecessors(nd)), Roadm)]
+    preamp_nodes_nobuild_uid = [nd.uid for nd in network.nodes() \
+        if isinstance(nd, Edfa) and isinstance(next(network.successors(nd)), Roadm)]
+    amp_nodes_nobuild = {nd.uid : nd for nd in network.nodes() \
+        if isinstance(nd, Edfa) and isinstance(next(network.predecessors(nd)), Roadm)}
+    preamp_nodes_nobuild = {nd.uid : nd for nd in network.nodes() \
+        if isinstance(nd, Edfa) and isinstance(next(network.successors(nd)), Roadm)}
+    # roadm dict with restrictions before build
+    roadms = {nd.uid: nd for nd in network.nodes() if isinstance(nd, Roadm)}
+    # Build the network once using the default power defined in SI in eqpt config
+    # power density : db2linp(ower_dbm": 0)/power_dbm": 0 * nb channels as defined by
+    # spacing, f_min and f_max
+    p_db = equipment['SI']['default'].power_dbm
+    p_total_db = p_db + lin2db(automatic_nch(equipment['SI']['default'].f_min,\
+        equipment['SI']['default'].f_max, equipment['SI']['default'].spacing))
+
+    build_network(network, equipment, p_db, p_total_db)
+
+    amp_nodes = [nd for nd in network.nodes() \
+        if isinstance(nd, Edfa) and isinstance(next(network.predecessors(nd)), Roadm)\
+           and next(network.predecessors(nd)).restrictions['booster_variety_list']]
+
+    preamp_nodes = [nd for nd in network.nodes() \
+        if isinstance(nd, Edfa) and isinstance(next(network.successors(nd)), Roadm)\
+           and next(network.successors(nd)).restrictions['preamp_variety_list']]
+
+    # check that previously existing amp are not changed
+    for amp in amp_nodes:
+        if amp.uid in amp_nodes_nobuild_uid:
+            print(amp.uid, amp.params.type_variety)
+            if not amp.params.type_variety == amp_nodes_nobuild[amp.uid].params.type_variety:
+                raise AssertionError()
+    for amp in preamp_nodes:
+        if amp.uid in preamp_nodes_nobuild_uid:
+            if not amp.params.type_variety == preamp_nodes_nobuild[amp.uid].params.type_variety:
+                raise AssertionError()
+    # check that restrictions are correctly applied
+    for amp in amp_nodes:
+        if amp.uid not in amp_nodes_nobuild_uid:
+            # and if roadm had no restrictions before build:
+            if restrictions['booster_variety_list'] and \
+               not roadms[next(network.predecessors(amp)).uid]\
+                         .restrictions['booster_variety_list']:
+                if not amp.params.type_variety in restrictions['booster_variety_list']:
+
+                    raise AssertionError()
+    for amp in preamp_nodes:
+        if amp.uid not in preamp_nodes_nobuild_uid:
+            if restrictions['preamp_variety_list'] and\
+            not roadms[next(network.successors(amp)).uid].restrictions['preamp_variety_list']:
+                if not amp.params.type_variety in restrictions['preamp_variety_list']:
+                    raise AssertionError()