Add script to auto-generate theoretical formula

Python/Theory_N_receivers_auto_develop.py

+import math
+import matplotlib.pyplot as plt
+import numpy as np
+import itertools
+
+VERBOSE = False
+
+
+def _print(*args):
+    if VERBOSE:
+        print(*args)
+
+
+def sgn(s, i):
+    """
+    Return 1 if i is in the list s, returns -1 otherwise
+    """
+    return 2 * (i in s) - 1
+
+
+def inter(L1, L2):
+    """
+    Return the intersection between two lists
+    """
+    return set(L1).intersection(L2)
+
+
+def Q(x):
+    """
+    Implementation of the Q-function
+    """
+    return (1 - math.erf(x / math.sqrt(2))) / 2
+
+
+def boundary_string(n, signs, epsilon, N_u):
+    """
+    Auto-computes the boundaries for the noise considering a user n 
+    """
+    # Compute, for each user, the coefficient in front of the power
+    # This can help simplify the equations a bit
+    power_coeffs = [0] * N_u
+
+    for i in range(1, n + 1):
+        if i in signs:
+            power_coeffs[i - 1] -= 1  # - \\sqrt{\\frac{P_{" + str(i) + "}}{2}}
+        else:
+            power_coeffs[i - 1] += 1  # + \\sqrt{\\frac{P_{" + str(i) + "}}{2}}
+
+    for i in epsilon:
+        if i in signs:
+            power_coeffs[i - 1] -= 2  # - 2 \\sqrt{\\frac{P_{" + str(i) + "}}{2}}
+        else:
+            power_coeffs[i - 1] += 2  # + 2 \\sqrt{\\frac{P_{" + str(i) + "}}{2}}
+
+    # Build the formula based on the contents of power_coeffs
+    formula = ""
+    nbr_items = 0
+    for i in range(N_u, 0, -1):
+        nbr_items += power_coeffs[i - 1] != 0
+
+        if power_coeffs[i - 1] == 2:
+            formula += " + 2 \\sqrt{\\frac{P_{" + str(i) + "}}{2}}"
+        elif power_coeffs[i - 1] == 1:
+            formula += " + \\sqrt{\\frac{P_{" + str(i) + "}}{2}}"
+        elif power_coeffs[i - 1] == -1:
+            formula += " - \\sqrt{\\frac{P_{" + str(i) + "}}{2}}"
+        elif power_coeffs[i - 1] == -2:
+            formula += " - 2 \\sqrt{\\frac{P_{" + str(i) + "}}{2}}"
+
+    # Format it (remove the leading " + ", and replace the leading " - " by "-")
+    if len(formula) > 0:
+        if formula[1] == "+":
+            formula = formula[3:]
+        else:
+            formula = "-" + formula[3:]
+
+    # Add parenthesis if necessary
+    if nbr_items > 2:
+        formula = "\\left(" + formula + "\\right)"
+
+    return formula
+
+
+def theory(n, g, sigma, P, N_u):
+    """
+    Auto-develops the analytical formula for user n, and computes its value under the given conditions
+    """
+    # All the powers must be > 0
+    if any(x < 0 for x in P):
+        raise ValueError("Every user's power must be > 0")
+
+    # Sigma must be > 0
+    if sigma <= 0:
+        raise ValueError("The noise's variance must be > 0")
+
+    # There must be at least 1 user
+    if N_u < 1:
+        raise ValueError("There must be at least 1 user")
+
+    if len(P) != N_u:
+        raise ValueError("Number of users and power values are incoherent")
+
+    # Attenuation cannot be <= 0
+    if g <= 0:
+        raise ValueError("Attenuation must be > 0")
+
+    # Generate all the possibilities for the signs
+    S = list(itertools.product([0, 1], repeat=N_u))
+
+    # Generatie all the possibilities of errors for n
+    E_n = list(itertools.product([0, 1], repeat=(N_u - n)))
+
+    # List of probabilities under each condition
+    Probas = []
+
+    # Variable storing the analytical formula
+    formula_string = ""  # Latex formula with P(...)
+    formula_string_Q = ""  # Latex formula with Q(...)
+    numbered_formula_string = ""
+    numbers_string = ""
+
+    for s in S:
+        # Find the indices with a positive sign
+        signs = [i + 1 for i in range(len(s)) if s[i] == 1]
+
+        # Find the indices with a negative sign
+        signs_bar = [i + 1 for i in range(len(s)) if s[i] == 0]
+
+        # Knowing the signs, the r_n can be computer
+        r = [math.sqrt(P[i] / 2) * sgn(signs, i + 1) for i in range(len(P))]
+
+        for e in E_n:
+            # Find the indices where an error is expected
+            epsilon = [n] + [n + 1 + i for i in range(len(e)) if e[i] == 1]
+
+            # Find the indices where there is no error
+            epsilon_bar = [n + 1 + i for i in range(len(e)) if e[i] == 0]
+
+            _print(f"-- S={signs}, !S={signs_bar}, ε={epsilon}, !ε={epsilon_bar} --")
+
+            # List of constraints for the upper boundary of Re(n/g)
+            Csup = []
+            Csup_formulas = []
+
+            # List of constraints for the lower boundary of Re(n/g)
+            Cinf = []
+            Cinf_formulas = []
+
+            # Compute the boundaries
+            _print(f"S∩ε={inter(epsilon, signs)}")
+            for i in inter(epsilon, signs):
+                epsilon_i = inter(epsilon, range(i + 1, N_u + 1))
+                Csup.append(-sum(r[:i]) - 2 * sum([r[j - 1] for j in epsilon_i]))
+
+                # Build the string describing this boundary
+                Csup_formulas.append(boundary_string(i, signs, epsilon_i, N_u))
+
+            _print(f"S∩!ε={inter(epsilon_bar, signs)}")
+            for i in inter(epsilon_bar, signs):
+                epsilon_i = inter(epsilon, range(i + 1, N_u + 1))
+                Cinf.append(-sum(r[:i]) - 2 * sum([r[j - 1] for j in epsilon_i]))
+
+                # Build the string describing this boundary
+                Cinf_formulas.append(boundary_string(i, signs, epsilon_i, N_u))
+
+            _print(f"!S∩ε={inter(epsilon, signs_bar)}")
+            for i in inter(epsilon, signs_bar):
+                epsilon_i = inter(epsilon, range(i + 1, N_u + 1))
+                Cinf.append(-sum(r[:i]) - 2 * sum([r[j - 1] for j in epsilon_i]))
+
+                # Build the string describing this boundary
+                Cinf_formulas.append(boundary_string(i, signs, epsilon_i, N_u))
+
+            _print(f"!S∩!ε={inter(epsilon_bar, signs_bar)}")
+            for i in inter(epsilon_bar, signs_bar):
+                epsilon_i = inter(epsilon, range(i + 1, N_u + 1))
+                Csup.append(-sum(r[:i]) - 2 * sum([r[j - 1] for j in epsilon_i]))
+
+                # Build the string describing this boundary
+                Csup_formulas.append(boundary_string(i, signs, epsilon_i, N_u))
+
+            _print(Cinf, Csup)
+            # Compute the value of the Q-function base on these contraints
+            if len(Cinf) == 0 and len(Csup) == 0:
+                formula_string += " + 1"
+                formula_string_Q += " + 1"
+                numbered_formula_string += " + 1"
+                Probas.append(1 / 2**N_u)
+            elif len(Cinf) == 0:
+                min_sup = min(Csup)
+                index_min_sup = Csup.index(min_sup)
+
+                formula_string += " + P\\left(\\frac{\\eta}{g} < " + Csup_formulas[index_min_sup] + "\\right)"
+                formula_string_Q += " + 1 - Q\\left(\\frac{g}{\\sigma} \\left(" +  Csup_formulas[index_min_sup] + "\\right)\\right)"
+                numbered_formula_string += f" + ℙ(η/g < {min_sup})"
+                Probas.append(1 / 2**N_u * (1 - Q(min_sup * g / sigma)))
+            elif len(Csup) == 0:
+                max_inf = max(Cinf)
+                index_max_inf = Cinf.index(max_inf)
+
+                formula_string += f" + P\\left(" + Cinf_formulas[index_max_inf] + " < \\frac{\\eta}{g}\\right)"
+                formula_string_Q += " + Q\\left(\\frac{g}{\\sigma} \\left(" +  Cinf_formulas[index_max_inf] + "\\right)\\right)"
+                numbered_formula_string += f" + ℙ({max_inf} < η/g)"
+                Probas.append(1 / 2**N_u * Q(max_inf * g / sigma))
+            else:
+                max_inf = max(Cinf)
+                min_sup = min(Csup)
+                index_max_inf = Cinf.index(max_inf)
+                index_min_sup = Csup.index(min_sup)
+
+                formula_string += f" + P\\left(" + Cinf_formulas[index_max_inf] + " < \\frac{\\eta}{g} < " + Csup_formulas[index_min_sup] + "\\right)"
+                numbered_formula_string += f" + ℙ({max_inf} < η/g < {min_sup})"
+
+                if max_inf <= min_sup:
+                    Probas.append(1 / 2**N_u * (Q(max_inf * g / sigma) - Q(min_sup * g / sigma)))
+                    formula_string_Q += " + Q\\left(\\frac{g}{\\sigma} \\left(" + Cinf_formulas[index_max_inf] + "\\right) \\right) - Q\\left(\\frac{g}{\\sigma} \\left(" + Csup_formulas[index_min_sup] + "\\right) \\right)"
+                else:
+                    Probas.append(0)
+                    formula_string_Q += " + 0"
+
+            numbers_string += f" + {Probas[-1]}"
+
+    formula_string = formula_string[3:]
+    formula_string = "\\frac{1}{" + str(2**N_u) + "} \\left(" + formula_string + "\\right)"
+
+    formula_string_Q = formula_string_Q[3:]
+    formula_string_Q = "\\frac{1}{" + str(2**N_u) + "} \\left(" + formula_string_Q + "\\right)"
+
+    numbered_formula_string = numbered_formula_string[3:]
+    numbered_formula_string = f"1/{2**N_u} ({numbered_formula_string})"
+
+    numbers_string = numbers_string[3:]
+
+    return formula_string, formula_string_Q, numbered_formula_string, numbers_string, sum(Probas)
+
+
+if __name__ == '__main__':
+    N_u = 2
+    n = 2
+    sigma = 1
+    g = 1
+    P = [10**k for k in range(N_u)]
+    print("N_u =", N_u)
+    print("Powers =", P)
+    print("Attenuation g =", g)
+    print("Noise with a variance of σ² =", sigma**2)
+
+    latex_formula, latex_formula_Q, formula_with_numbers, sum_of_numbers, result = theory(n, g, sigma, P, N_u)
+    print(latex_formula)
+    print("=", latex_formula_Q)
+    print("=", formula_with_numbers)
+    print("≈", sum_of_numbers)
+    print("≈", result)