Mathematica implementation of Mason's gain formula

Mason's gain formula (MGF) is a method for finding the transfer function of a linear signal-flow graph (SFG).

I'm trying to port the MATLAB code to Mathematica code, but the Mathematica code outputs false result.

(* ::Package:: *)

(* masonFormula.m
   Translates the mason.m Matlab function to Mathematica.
   Calculates the transfer function (Numerator/Denominator) for a
   signal flow graph described in a .net file, using Mason's Gain Formula.

   Based on:
   mason.m by Rob Walton
   TRLabs and the University of Calgary, Alberta, Canada
   Date: January 25th 1999
   Revised: January 20th 2000

   USAGE:
     {Numerator, Denominator} = masonFormula["Netfile", StartNode, StopNode]

     Netfile     - is a string with the name of the netfile to load
     StartNode   - is the integer number describing the independent input node
     StopNode    - is the integer number describing the dependent output node
     Numerator   - is a symbolic expression for the Numerator
     Denominator - is a symbolic expression for the Denominator
*)


(* --- Helper Function Definitions --- *)

(* Function to check if two lists of nodes have any common elements *)
AreTouchingLoops[nodes1_List, nodes2_List] :=
  Length[Intersection[nodes1, nodes2]] > 0;

(* Function to check if two lists of coefficients represent the same loop OR
   the same combination of non-touching loops.
   This means they contain the same multiset of non-zero coefficient numbers *)
IsSameLoopCombination[coeffs1_List, coeffs2_List] :=
  Sort[Select[coeffs1, # > 0 &]] == Sort[Select[coeffs2, # > 0 &]];

(* Recursive function to find simple paths (no repeated nodes) *)
(* Returns a list of {pathCoefficients, pathNodes} pairs *)
findPaths[curr_, stop_, pathCoeffs_, visitedNodes_, net_] :=
  Module[{newVisitedNodes, outgoingEdges, results},
    (* Check if current node has been visited in this path already *)
    If[MemberQ[visitedNodes, curr], Return[{}]];

    newVisitedNodes = Append[visitedNodes, curr];

    (* Base case: Reached the stop node *)
    If[curr == stop,
      (* A path must have at least one edge *)
      If[Length[pathCoeffs] > 0,
        Return[{{pathCoeffs, newVisitedNodes}}],
        Return[{}] (* Start == Stop with no edges, not a valid path/loop in this context *)
      ]
    ];

    (* Find outgoing edges from the current node *)
    outgoingEdges = Select[net, #[[2]] == curr &]; (* Edges where start node is curr *)

    (* Recursively call for each outgoing edge *)
    results = Flatten[
      Table[
        (* The edge format is {coeffIndex, startNode, endNode} *)
        findPaths[edge[[3]], stop, Append[pathCoeffs, edge[[1]]], newVisitedNodes, net],
        {edge, outgoingEdges}
      ],
      1 (* Flatten by one level *)
    ];

    Return[results];
  ];

(* Function to remove duplicate loop/combination descriptions from a list *)
RemoveDuplicateLoopCombinations[loopList_List] :=
  Module[{uniqueLoops = {}},
    Scan[
      (* Use IsSameLoopCombination for comparison *)
      If[! MemberQ[uniqueLoops, #, (IsSameLoopCombination[#1["Coeff"], #2["Coeff"]] &)],
         AppendTo[uniqueLoops, #]
      ] &,
      loopList
    ];
    Return[uniqueLoops];
  ];

(* Function to convert a list of coefficient numbers into a symbolic product *)
CoeffsToProduct[coefficients_List, coeffSymbols_List] :=
  Module[{validCoeffs},
    validCoeffs = Select[coefficients, # > 0 &]; (* Remove padding/zeros *)
    If[Length[validCoeffs] == 0,
      1 (* Product of empty set is 1. Represents the '1' term in Delta and Delta_k *)
    ,
      (* Map coefficient numbers (indices) to symbols and take the product *)
      Apply[Times, coeffSymbols[[validCoeffs]]]
    ]
  ];

(* Function to calculate the sum of products for loop combinations of a given order
   that do not touch a given set of nodes (pNodes).
   'L' is the Association holding loop combinations grouped by order.
   'order' specifies which order combinations to consider (e.g., 1 for single loops, 2 for pairs, etc.).
   'pNodes' is the list of nodes that the loop combination should NOT touch.
   Returns the sum: Sum[ Product[coeffs] ] for all valid non-touching combinations of the given order.
*)
SumOfProductsNonTouching[L_Association, order_Integer, pNodes_List, coeffSymbols_List] :=
  Module[{orderLoopCombinations, nonTouchingCombinations, products},
    (* Check if combinations of this order exist *)
    If[!KeyExistsQ[L, order] || Length[L[order]] == 0,
      Return[0] (* No combinations of this order *)
    ];

    orderLoopCombinations = L[order];

    (* Find loop combinations of this order that do not touch the path nodes *)
    (* #["Node"] contains ALL nodes from the constituent loops in the combination *)
    nonTouchingCombinations = Select[orderLoopCombinations, ! AreTouchingLoops[#["Node"], pNodes] &];

    (* Calculate the product of coefficients for each non-touching combination *)
    (* #["Coeff"] contains ALL coefficients from the constituent loops *)
    products = CoeffsToProduct[#["Coeff"], coeffSymbols] & /@ nonTouchingCombinations;

    (* Sum the products *)
    Total[products]
  ];


(* --- Main Function --- *)

masonFormula[netFile_String, start_Integer, stop_Integer] :=
  Module[{
    file, line, coeff, net = {}, coeffNames = {}, numCoeffs, numNodes,
    rawPaths, paths = {},
    rawLoops = {}, loopCombinations = <||>, (* Renamed 'loops' to 'loopCombinations' for clarity *)
    coeffSymbols,
    numeratorExpr, denominatorExpr,
    allNodes
    },

    (* --- Load the net description file --- *)
    file = OpenRead[netFile];
    If[file === $Failed,
  Print["*** File, ", netFile, ", not found ***"];
  Return[$Failed]
    ];

    While[True,
      line = ReadList[file, Number, 3];
      If[line === EndOfFile, Break[]];
      coeff = Read[file, Word];
      If[coeff === EndOfFile, Break[]];
      AppendTo[net, line];
      AppendTo[coeffNames, coeff];
      ReadList[file, __, RecordSeparators -> {"\n"}];
    ];
    Close[file];

    numCoeffs = Length[net];
    If[numCoeffs === 0,
      Print["*** No data read from file ", netFile, " ***"];
      Return[$Failed]
    ];

    (* Determine the set of all nodes involved *)
    allNodes = Union[Flatten[net[[All, {2, 3}]]]];
    numNodes = Length[allNodes]; (* Max node number might be different *)

    (* Create symbolic variables from coefficient names *)
    coeffSymbols = Symbol /@ coeffNames;

    (* --- Find all forward paths connecting start and stop nodes --- *)
    rawPaths = findPaths[start, stop, {}, {}, net];

    If[Length[rawPaths] == 0,
      Print["*** There are no paths connecting nodes ", start, " and ", stop, " ***"];
      (* According to Mason's formula, if there are no paths, the gain is 0.
         The denominator (Delta) is usually non-zero if loops exist, or 1 otherwise. *)
      (* We still need to calculate Delta for completeness, or just return 0/1 *)
       (* Let's proceed to calculate Delta anyway, and return {0, Delta} *)
       (* numeratorExpr = 0; (* Skip path processing *) *)
    ];

    (* Structure path data *)
    paths = Map[<|"Coeff" -> #[[1]], "Node" -> #[[2]]|> &, rawPaths];


    (* --- Determine all first-Order loops (L_1 term) --- *)
    (* Find loops starting from *any* node involved in the network *)
    For[i = 1, i <= Length[allNodes], ++i,
        node = allNodes[[i]];
        rawLoops = Join[rawLoops, findPaths[node, node, {}, {}, net]];
    ];


    (* Remove duplicate loops using coefficient comparison *)
    (* findPaths can return the same loop starting from different nodes *)
    rawLoops = RemoveDuplicateLoopCombinations[rawLoops]; (* Use the renamed function *)

    (* Structure loop data for order 1 *)
    If[Length[rawLoops] > 0,
        loopCombinations[1] = Map[<|"Coeff" -> #[[1]], "Node" -> #[[2]]|> &, rawLoops];
    ,
        loopCombinations[1] = {}; (* Ensure key exists even if no loops *)
    ];


    (* --- Determine higher-order loop combination terms (L_m for m > 1) --- *)
    (* These represent combinations of non-touching loops *)
    order = 1;
    While[True,
      order++;
      loopCombinations[order] = {}; (* Initialize list for this order *)

      (* Need loops of order 1 and combinations of order (order-1) to build order 'order' *)
      If[!KeyExistsQ[loopCombinations, 1] || Length[loopCombinations[1]] == 0 ||
         !KeyExistsQ[loopCombinations, order - 1] || Length[loopCombinations[order - 1]] == 0,
         Break[] (* No loops of previous orders, so no higher order combinations *)
      ];

      (* Compare each first-order loop with each (order-1)-th combination *)
      For[i = 1, i <= Length[loopCombinations[1]], ++i,
        For[j = 1, j <= Length[loopCombinations[order - 1]], ++j,
          (* Check if the first-order loop touches the existing combination *)
          If[! AreTouchingLoops[loopCombinations[1, i, "Node"], loopCombinations[order - 1, j, "Node"]],
            (* Non-touching found, combine them *)
            combinedCoeffs = Join[loopCombinations[1, i, "Coeff"], loopCombinations[order - 1, j, "Coeff"]];
            combinedNodes = Join[loopCombinations[1, i, "Node"], loopCombinations[order - 1, j, "Node"]];

            (* Check if this exact combination (based on coeffs) is already in loopCombinations[order] *)
            isDuplicate = MemberQ[loopCombinations[order], <|"Coeff" -> combinedCoeffs|>,
                                 (IsSameLoopCombination[#1["Coeff"], #2["Coeff"]] &)];

            If[! isDuplicate,
              AppendTo[loopCombinations[order], <|"Coeff" -> combinedCoeffs, "Node" -> combinedNodes|>];
            ];
          ]
        ]
      ];

      (* If no combinations of this order were found, stop *)
      If[Length[loopCombinations[order]] == 0,
          KeyDropFrom[loopCombinations, order]; (* Remove empty entry *)
          Break[]
      ];
      (* Important: Remove duplicates *after* generating all possibilities for the order *)
      loopCombinations[order] = RemoveDuplicateLoopCombinations[loopCombinations[order]];
       If[Length[loopCombinations[order]] == 0, Break[]]; (* Check again after deduplication *)
    ];

    (* --- Display Info (Optional) --- *)
    Print["\n-- Network Info --"];
    Print["Net File   : ", netFile];
    Print["Start Node : ", start];
    Print["Stop Node  : ", stop];

    Print["\n----- Paths (P_k) -----"];
    If[Length[paths] == 0,
        Print["(No paths found)"],
        For[i = 1, i <= Length[paths], ++i,
             (* Print symbolic path gain *)
             Print["P", i, " (Nodes: ", paths[[i, "Node"]], ") Gain: ", CoeffsToProduct[paths[[i, "Coeff"]], coeffSymbols]];
        ]
    ];

    Print["\n----- Loop Combinations (Terms for Delta) -----"];
    If[Length[Keys[loopCombinations]] == 0 || Total[Length /@ Values[loopCombinations]] == 0,
        Print["(No loops found)"],
        Scan[
            Function[{orderKey},
                Print["\n- Order ", orderKey, " Combinations (Term L",orderKey,") -"];
                If[Length[loopCombinations[orderKey]] == 0,
                    Print["(None)"],
                    For[i = 1, i <= Length[loopCombinations[orderKey]], ++i,
                         (* Print symbolic product for this combination *)
                         Print["Term ", i, " (Coeffs: ", loopCombinations[orderKey, i, "Coeff"], ") Product: ", CoeffsToProduct[loopCombinations[orderKey, i, "Coeff"] , coeffSymbols] ];
                    ]
                ]
            ],
            Sort[Keys[loopCombinations]] (* Print orders in numerical sequence *)
        ]
    ];
    Print[""]; (* Add a newline *)


    (* --- Generate the final equations (Consistent with Mason's Gain Formula) --- *)

    (* Determine Denominator (Delta) *)
    (* Delta = 1 - L1 + L2 - L3 + ... *)
    denominatorExpr = 1;
    Scan[
      Function[{orderKey},
        sign = If[OddQ[orderKey], -1, 1];
        (* SumOfProductsNonTouching with empty node list gives sum of products for ALL combinations of this order (L_m) *)
        denominatorExpr += sign * SumOfProductsNonTouching[loopCombinations, orderKey, {}, coeffSymbols];
      ],
      Sort[Keys[loopCombinations]]
    ];

    (* Determine Numerator (Sum_k P_k * Delta_k) *)
    numeratorExpr = 0;
    If[Length[paths] > 0, (* Only calculate if paths exist *)
        For[pathn = 1, pathn <= Length[paths], ++pathn,
          pathData = paths[[pathn]];
          pathGain = CoeffsToProduct[pathData["Coeff"], coeffSymbols]; (* P_k *)
          pathNodes = pathData["Node"];

          (* Calculate Delta_k for this path *)
          (* Delta_k = 1 - L1(k) + L2(k) - L3(k) + ... *)
          deltaK = 1;
          Scan[
            Function[{orderKey},
              sign = If[OddQ[orderKey], -1, 1];
              (* SumOfProductsNonTouching with pathNodes gives sum of products for combinations NOT touching the path (L_m(k)) *)
              deltaK += sign * SumOfProductsNonTouching[loopCombinations, orderKey, pathNodes, coeffSymbols];
            ],
            Sort[Keys[loopCombinations]]
          ];

          numeratorExpr += pathGain * deltaK; (* Sum P_k * Delta_k *)
        ];
    ]; (* End If paths exist *)

    Print["The variables returned are symbolic expressions for the"];
    Print["Numerator and Denominator of the transfer function G = Numerator / Denominator."];
    Print["Use Simplify[Numerator/Denominator] to potentially simplify the result.\n"];

    Return[{numeratorExpr, denominatorExpr}];
  ];


(* --- Example Usage --- *)

(*
   1    1   2   s21
   2    4   2   s22
   3    4   3   s12
   4    1   3   s11
   5    2   4   R2
*)
(* Create the example.net file *)
Export["example.net",
"1  1   2   s21
2   4   2   s22
3   4   3   s12
4   1   3   s11
5   2   4   R2", "Text"];


(* Calculate transfer function from node 1 to node 3 *)
 {num, den} = masonFormula["example.net", 1, 3]; 

(* Print the results *)
 Print["Numerator: ", num]; 
 Print["Denominator: ", den]; 

(* Print simplified transfer function *)
 Print["Transfer Function (Simplified): ", Simplify[num/den]]; 

My question is: how to fix the wrong code? Now the finding paths part is correct, whereas the finding loop part is not correct. Thanks in advance.