<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Optimization of Degree Distributions on the BEC
This code is provided as supplementary material of the lecture Channel Coding 2 - Advanced Methods.
This code illustrates
* Using linear programming to optimize degree distributions on the BEC
Step1:  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}} 
Step2:  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.
Step3:  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$
Step4:  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.
Step5:  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.
Step6:  Run binary search to find best irregular code for a given target rate on the BEC.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plot
from ipywidgets import interactive
import ipywidgets as widgets
import math
from pulp import *
%matplotlib inline  
Explanation: Optimization of Degree Distributions on the BEC
This code is provided as supplementary material of the lecture Channel Coding 2 - Advanced Methods.
This code illustrates
* Using linear programming to optimize degree distributions on the BEC
End of explanation
# 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
Explanation: 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.
End of explanation
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
Explanation: 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.
End of explanation
best_lambda = find_best_lambda(0.2949219, 16, 12.98)
print(np.poly1d(best_lambda, variable='Z'))
Explanation: 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$
End of explanation
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
Explanation: 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.
End of explanation
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
Explanation: 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.
End of explanation
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
Explanation: Run binary search to find best irregular code for a given target rate on the BEC.
End of explanation
<END_TASK>