In [1]:
import numpy as np
import matplotlib.pyplot as plot
from ipywidgets import interactive
import ipywidgets as widgets
import math
from pulp import *
%matplotlib inline
We specify the check node degree distribution polynomial $\rho(Z)$ by fixing the average check node degree $d_{\mathtt{c},\text{avg}}$ and assuming that the code contains only check nodes with degrees $\tilde{d}_{\mathtt{c}} := \lfloor d_{\mathtt{c},\text{avg}}\rfloor$ and $\tilde{d}_{\mathtt{c}}+1$. This is the so-called check-concentrated degree distribution. As shown in the lecture, we have: $$ \rho(Z) = \frac{\tilde{d}_{\mathtt{c}}(\tilde{d}_{\mathtt{c}}+1-d_{\mathtt{c},\text{avg}})}{d_{\mathtt{c},\text{avg}}}Z^{\tilde{d}_{\mathtt{c}}-1} + \frac{d_{\mathtt{c},\text{avg}}-\tilde{d}_{\mathtt{c}}(\tilde{d}_{\mathtt{c}}+1-d_{\mathtt{c},\text{avg}})}{d_{\mathtt{c},\text{avg}}}Z^{\tilde{d}_{\mathtt{c}}} $$ The following function converts $d_{\mathtt{c},\text{avg}}$ into a polynomial $\rho(Z)$ which is given as an array where the first entry corresponds to the largest exponents and the last entry corresponds to the constant part.
In [2]:
# returns rho polynomial (highest exponents first) corresponding to average check node degree c_avg
def c_avg_to_rho(c_avg):
ct = math.floor(c_avg)
r1 = ct*(ct+1-c_avg)/c_avg
r2 = (c_avg - ct*(ct+1-c_avg))/c_avg
rho_poly = np.concatenate(([r2,r1], np.zeros(ct-1)))
return rho_poly
The following function solves the optimization problem that returns the best $\lambda(Z)$ for a given BEC erasure probability $\epsilon$, for an average check node degree $d_{\mathtt{c},\text{avg}}$, and for a maximum variable node degree $d_{\mathtt{v},\max}$. This optimization problem is derived in the lecture as $$ \begin{aligned} & \underset{\lambda_1,\ldots,\lambda_{d_{\mathtt{v},\max}}}{\text{maximize}} & & \sum_{i=1}^{d_{\mathtt{v},\max}}\frac{\lambda_i}{i} \\ & \text{subject to} & & \lambda_1 = 0 \\ & & & \lambda_i \geq 0, \quad \forall i \in\{2,3,\ldots,d_{\mathtt{v},\max}\} \\ & & & \sum_{i=2}^{d_{\mathtt{v},\max}}\lambda_i = 1 \\ & & & \sum_{i=2}^{d_{\mathtt{v},\max}}\lambda_i\cdot \epsilon(1-\rho(1-\tilde{\xi}_j))^{i-1}-\tilde{\xi}_j \leq 0,\quad \forall j \in \{1,\ldots, D\} \\ & & & \lambda_2 \leq \frac{1}{\epsilon\rho^\prime(1)} = \frac{1}{\epsilon\sum_{i=2}^{d_{\mathtt{c},\max}}(i-1)\rho_i} \end{aligned} $$
If this optimization problem is feasible, then the function returns the polynomial $\lambda(Z)$ as a coefficient array where the first entry corresponds to the largest exponent ($\lambda_{d_{\mathtt{v},\max}}$) and the last entry to the lowest exponent ($\lambda_1$). If the optimization problem has no solution (e.g., it is unfeasible), then the empty vector is returned.
In [3]:
def find_best_lambda(epsilon, v_max, c_avg):
rho = c_avg_to_rho(c_avg)
# quantization of EXIT chart
D = 500
xi_range = np.arange(1.0, D+1, 1)/D
# Linear Programming model, maximize target expression
model = pulp.LpProblem("Finding best lambda problem", pulp.LpMaximize)
# definition of variables, v_max entries \lambda_i that are between 0 and 1 (implicit declaration of constraint 2)
v_lambda = pulp.LpVariable.dicts("lambda", range(v_max),0,1)
# objective function
cv = 1/np.arange(v_max,0,-1)
model += pulp.lpSum(v_lambda[i]*cv[i] for i in range(v_max))
# constraints
# constraint 1, no variable nodes of degree 1
model += v_lambda[v_max-1] == 0
# constraint 3, sum of lambda_i must be 1
model += pulp.lpSum(v_lambda[i] for i in range(v_max))==1
# constraints 4, fixed point condition for all the descrete xi values (a total number of D, for each \xi)
for xi in xi_range:
model += pulp.lpSum(v_lambda[j] * epsilon * (1-np.polyval(rho,1.0-xi))**(v_max-1-j) for j in range(v_max))-xi <= 0
# constraint 5, stability condition
model += v_lambda[v_max-2] <= 1/epsilon/np.polyval(np.polyder(rho),1.0)
model.solve()
if model.status != 1:
r_lambda = []
else:
r_lambda = [v_lambda[i].varValue for i in range(v_max)]
return r_lambda
As an example, we consider the case of optimization carried out in the lecture after 9 iterations, where we have $\epsilon = 0.2949219$ and $d_{\mathtt{c},\text{avg}} = 12.98$ with $d_{\mathtt{v},\max}=16$
In [4]:
best_lambda = find_best_lambda(0.2949219, 16, 12.98)
print(np.poly1d(best_lambda, variable='Z'))
In the following, we provide an interactive widget that allows you to choose the parameters of the optimization yourself and get the best possible $\lambda(Z)$. Additionally, the EXIT chart is plotted to visualize the good fit of the obtained degree distribution.
In [8]:
def best_lambda_interactive(epsilon, c_avg, v_max):
# get lambda and rho polynomial from optimization and from c_avg, respectively
p_lambda = find_best_lambda(epsilon, v_max, c_avg)
p_rho = c_avg_to_rho(c_avg)
# if optimization successful, compute rate and show plot
if not p_lambda:
print('Optimization infeasible, no solution found')
else:
design_rate = 1 - np.polyval(np.polyint(p_rho),1)/np.polyval(np.polyint(p_lambda),1)
if design_rate <= 0:
print('Optimization feasible, but no code with positive rate found')
else:
print("Lambda polynomial:")
print(np.poly1d(p_lambda, variable='Z'))
print("Design rate r_d = %1.3f" % design_rate)
# Plot EXIT-Chart
print("EXIT Chart:")
plot.figure(3)
x = np.linspace(0, 1, num=100)
y_v = [1 - epsilon*np.polyval(p_lambda, 1-xv) for xv in x]
y_c = [np.polyval(p_rho,xv) for xv in x]
plot.plot(x, y_v, '#7030A0')
plot.plot(y_c, x, '#008000')
plot.axis('equal')
plot.gca().set_aspect('equal', adjustable='box')
plot.xlim(0,1)
plot.ylim(0,1)
plot.xlabel('$I^{[A,V]}$, $I^{[E,C]}$')
plot.ylabel('$I^{[E,V]}$, $I^{[A,C]}$')
plot.grid()
plot.show()
interactive_plot = interactive(best_lambda_interactive, \
epsilon=widgets.FloatSlider(min=0.01,max=1,step=0.001,value=0.5, continuous_update=False, description=r'\(\epsilon\)',layout=widgets.Layout(width='50%')), \
c_avg = widgets.FloatSlider(min=3,max=20,step=0.1,value=4, continuous_update=False, description=r'\(d_{\mathtt{c},\text{avg}}\)'), \
v_max = widgets.IntSlider(min=3, max=20, step=1, value=16, continuous_update=False, description=r'\(d_{\mathtt{v},\max}\)'))
output = interactive_plot.children[-1]
output.layout.height = '400px'
interactive_plot
Now, we carry out the optimization over a wide range of $d_{\mathtt{c},\text{avg}}$ values for a given $\epsilon$ and find the largest possible rate.
In [9]:
def find_best_rate(epsilon, v_max, c_max):
c_range = np.linspace(3, c_max, num=100)
rates = np.zeros_like(c_range)
# loop over all c_avg, add progress bar
f = widgets.FloatProgress(min=0, max=np.size(c_range))
display(f)
for index,c_avg in enumerate(c_range):
f.value += 1
p_lambda = find_best_lambda(epsilon, v_max, c_avg)
p_rho = c_avg_to_rho(c_avg)
if p_lambda:
design_rate = 1 - np.polyval(np.polyint(p_rho),1)/np.polyval(np.polyint(p_lambda),1)
if design_rate >= 0:
rates[index] = design_rate
# find largest rate
largest_rate_index = np.argmax(rates)
best_lambda = find_best_lambda(epsilon, v_max, c_range[largest_rate_index])
print("Found best code of rate %1.3f for average check node degree of %1.2f" % (rates[largest_rate_index], c_range[largest_rate_index]))
print("Corresponding lambda polynomial")
print(np.poly1d(best_lambda, variable='Z'))
# Plot curve with all obtained results
plot.figure(4, figsize=(10,3))
plot.plot(c_range, rates, 'b')
plot.plot(c_range[largest_rate_index], rates[largest_rate_index], 'bs')
plot.xlim(3, c_max)
plot.ylim(0, (1.1*(1-epsilon)))
plot.xlabel('$d_{c,avg}$')
plot.ylabel('design rate $r_d$')
plot.grid()
plot.show()
return rates[largest_rate_index]
interactive_optim = interactive(find_best_rate, \
epsilon=widgets.FloatSlider(min=0.01,max=1,step=0.001,value=0.5, continuous_update=False, description=r'\(\epsilon\)',layout=widgets.Layout(width='50%')), \
v_max = widgets.IntSlider(min=3, max=20, step=1, value=16, continuous_update=False, description=r'\(d_{\mathtt{v},\max}\)'), \
c_max = widgets.IntSlider(min=3, max=40, step=1, value=22, continuous_update=False, description=r'\(d_{\mathtt{c},\max}\)'))
output = interactive_optim.children[-1]
output.layout.height = '400px'
interactive_optim
Run binary search to find best irregular code for a given target rate on the BEC.
In [9]:
target_rate = 0.7
dv_max = 16
dc_max = 22
T_Delta = 0.001
epsilon = 0.5
Delta_epsilon = 0.5
while Delta_epsilon >= T_Delta:
print('Running optimization for epsilon = %1.5f' % epsilon)
rate = find_best_rate(epsilon, dv_max, dc_max)
if rate > target_rate:
epsilon = epsilon + Delta_epsilon / 2
else:
epsilon = epsilon - Delta_epsilon / 2
Delta_epsilon = Delta_epsilon / 2
In [ ]: