# LDPC_simulation_stats_plots.m

Jan 23rd, 2022
1,143
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. function [H, num_variables, num_checks, num_edges, rate_expected, counter_tries, initial_errors, final_errors]= LDPC_simulation_stats_plots(n, l_optimum, r_optimum, PROCESS)
2.
3.     %% Estw oti brhkame ta PANTA SWSTA
4.     l_max = length(l_optimum);
5.     r_max = length(r_optimum);
6.
7.     % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
8.     % MPOREI NA MIN BGOUN AKERAIOI!!!
9.
10.     % l_optimum = [0.3, 0.5, 0.2];
11.     % r_optimum = [0.2, 0.2, 0.1, 0.1, 0.3, 0.1];
12.     % n = 200;
13.
14.
15.     %% Create Pi dia i and Lamda i dia i
16.     r_i_dia_i_list = zeros(1, r_max);
17.     l_i_dia_i_list = zeros(1, l_max);
18.     for index = 1:r_max
19.         r_i_dia_i_list(index) = r_optimum(index) / (index + 1);     % Epeidi arxizei apo to 2 enw emeis apo to 1
20.     end
21.
22.     for index = 1:l_max
23.         l_i_dia_i_list(index) = l_optimum(index) / (index + 1);
24.     end
25.     SUM = sum(l_i_dia_i_list);
26.
27.     r_i_dia_i_list;
28.     l_i_dia_i_list;
29.
30.     %% Create LAMDA AND RHO LISTS
31.     RHO = r_i_dia_i_list / SUM * n;
32.     LAMDA = l_i_dia_i_list / SUM * n;
33.     RHO_HAT = floor(RHO);
34.     LAMDA_HAT = floor(LAMDA);
35.     % Calculate A
36.     A_condition = sum(RHO - RHO_HAT);
37.
38.
39.     %% GENERAL CONSTRAINTS
40.     % Declare 2 sum constant variables
41.     constant_19 = n - sum(LAMDA_HAT);
42.     i_r = 2:r_max+1;
43.     i_l = 2:l_max+1;
44.     constant_20 = sum(i_r .* RHO_HAT) - sum(i_l .* LAMDA_HAT);
45.
46.     % Prwta stoixeia ta "l", kai meta ta "r"
47.     Aeq = zeros(2, l_max + r_max);
48.     Beq = zeros(2, 1);
49.
50.     % Constraint 19
51.     Aeq(1, 1:l_max) = 1;
52.     Beq(1, 1) = constant_19;
53.     % Constraint 20
54.     Aeq(2, 1:l_max) = i_l;
55.     Aeq(2, l_max+1 : l_max+r_max) = - i_r;
56.     Beq(2, 1) = constant_20;
57.     % Constraint 21
58.
59.
60.     %% CASE CONSTRAINTS
61.     % CASE 1 -  ANISOTIKOS - MATRIX A1
62.     A1 = zeros(1, l_max + r_max);
63.     A1(1, l_max+1 : l_max+r_max) = -1;
64.     B1 = -ceil(A_condition);
65.
66.     % CASE 2
67.     A2 = zeros(1, l_max + r_max);
68.     A2(1, l_max+1 : l_max+r_max) = 1;
69.     B2 = floor(A_condition);
70.
71.     % OBJECTIVE FUNCTIONS
72.     % CASE 1
73.     f1 = zeros(1, l_max + r_max);
74.     f1(1, l_max+1 : l_max+r_max) = 1;
75.     f1 = f1';
76.
77.     % CASE 2
78.     f2 = -f1;
79.
80.
81.     %% SIMULATIONS
82.     lb = zeros(1, l_max + r_max)';
83.     ub = ones(1, l_max + r_max)';
84.
85.     intcon = 1:l_max + r_max;
86.     % CASE 1
87.     result1 = intlinprog(f1, intcon, A1, B1 , Aeq, Beq, lb, ub);
88.     % CASE 2
89.     result2 = intlinprog(f2, intcon, A2, B2 , Aeq, Beq, lb, ub);
90.
91.     if isempty(result1) && isempty(result2)
92.         H = NaN;
93.         num_variables = NaN;
94.         num_checks = NaN;
95.         disp("Not possible solution found - exit code 1")
96.         return
97.     end
98. %     % CASE 1
99. %     result1 = linprog(f1, A1, B1 , Aeq, Beq, lb, ub);
100. %     % CASE 2
101. %     result2 = linprog(f2, A2, B2 , Aeq, Beq, lb, ub);
102.
103.     %% FIND THE BEST
104.
105.     if isempty(result1)
106.
107.         result = result2;
108.     elseif isempty(result2)
109.
110.         result = result1;
111.     else
112.          error1 = (sum(result1(l_max+1 : l_max+r_max)) - A_condition) ^ 2;
113.          error2 = (sum(result2(l_max+1 : l_max+r_max)) - A_condition) ^ 2;
114.          if error1 >= error2
115.             result = result2;
116.          else
117.             result = result1;
118.          end
119.     end
120.
121.
122.     % result = [xl1, xl2, xl3, xr1, ..., xr6] ---> length = 9
123.
124.
125.     %% FIND THE NEW LAMDA AND RHO - UPDATE
126.     result = floor(result');
127.     LAMDA_LIST = LAMDA_HAT + result(1 : l_max);
128.     RHO_LIST = RHO_HAT + result(l_max + 1 : l_max + r_max);
129.
130.
131.
132.
133.
134.
135.
136.
137.     %% DISPLAY AKMES
138.     % !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139.     % MPOREI KAPOIA NOUMERA NA MIN EINAI INTEGERS, ENW ANAPARISTOUN # AKMWN !
140.     num_variable_edges = sum(LAMDA_LIST .* i_l);
141.     num_check_edges = sum(RHO_LIST .* i_r);
142.     num_variables = sum(LAMDA_LIST);
143.     num_checks = sum(RHO_LIST);
144.     num_edges = num_variable_edges;
145.     n = num_variables;
146.     k = num_variables - num_checks;
147.     rate_expected = k / n;
148.
149.     % FOR CALLING THE FUNCTION
150.     counter_tries = 10^6;
151.     initial_errors = 10^7;
152.     final_errors = 10^8;
153.
154.
155.
156.     %   Above: L'(1) = P'(1)  (sxesh 2 kai 3)
157.
158.     if num_check_edges ~= num_variable_edges
159.         H = NaN;
160.         num_variables = NaN;
161.         num_checks = NaN;
162.         disp("Not possible solution found - exit code 2")
163.         return
164.     end
165.
166.
167.
168.
169.     % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
170.     % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
171.     % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
172.
173.     %% MATRIX 545 x 545 (about sockets)
174.     sockets = num_check_edges;
175.     checks = num_checks;
176.     vars = num_variables;
177.     comb1 = 1:sockets;
178.
179.     %% Create a vector with 82 2's and 91 3's and 27 4's (n=200)
180.     var_temp = zeros(1, vars);
181.     counter = 1;
182.     for i = 1:length(LAMDA_LIST)
183.         var_temp(counter : counter + LAMDA_LIST(i) - 1) = ...
184.         (i+1) * ones(1, LAMDA_LIST(i));
185.         counter = LAMDA_LIST(i) + counter;
186.     end
187.
188.     check_temp = zeros(1, checks);
189.     counter = 1;
190.     for i = 1:length(RHO_LIST)
191.         check_temp(counter : counter + RHO_LIST(i) - 1) = ...
192.         (i+1) * ones(1, RHO_LIST(i));
193.         counter = RHO_LIST(i) + counter;
194.     end
195.
196.     %% ACCUMULATIVE RANGES
197.     var_ranges = zeros(1, length(var_temp));
198.     var_ranges = [0 var_ranges];
199.
200.     for i = 1:length(var_temp)
201.         S = var_ranges(i) + var_temp(i);
202.         var_ranges(i+1) = S;
203.     end
204.
205.     check_ranges = zeros(1, length(check_temp));
206.     check_ranges = [0 check_ranges];
207.
208.     for i = 1:length(check_temp)
209.         S = check_ranges(i) + check_temp(i);
210.         check_ranges(i+1) = S;
211.     end
212.
213.
214.
215.
216.
217.
218.     anakyklwseis = 1000;
219.     H = zeros(checks, vars);
220.     counter_tries = 0;
221.
222.
223.
224.
225.
226.
227.
228.
229.     if PROCESS == "find_appropriate_n"
230.         %% FOR SIMULATIONS PURPOSE - DELETED LATER
231.         H = 1;
232.         num_variables = NaN;
233.         num_checks = NaN;
234.         disp("Possible solution found - exit code 3")
235.         return
236.
237.
238.
239.
240.
241.
242.
243.
244.     elseif PROCESS == "anakyklwseis"
245.
246.         display('*********************************************');
247.         display('*********************************************');
248.         pretty_display("n = ", n);
249.         pretty_display("Variable degrees = ", LAMDA_LIST);
250.         pretty_display("Check degrees = ", RHO_LIST);
251.         pretty_display("Num of variable edges = ", num_variable_edges);
252.         pretty_display("Num of check edges = ", num_check_edges);
253.         pretty_display("Num of VARIABLES = ", num_variables);
254.         pretty_display("Num of CHECKS = ", num_checks);
255.
256.         %% WHILE LOOP IN ORDER TO AVOID ANAKYKLWSEIS
257.         while anakyklwseis ~= 0
258.
259.             counter_tries = counter_tries + 1;
260.             socket_matrix = zeros(sockets);
261.             % Socket Matrix is 545 x 545
262.             comb2 = comb1(randperm(length(comb1)));
263.             for i = 1:sockets
264.                 ch = comb1(i);
265.                 var = comb2(i);
266.                 socket_matrix(ch, var) = 1;
267.             end
268.
269.
270.             %% Decomposition onto submatrices - Creation of matrix
271.             % var_ranges and check-ranges have 1 more element (0 in the beginning)
272.             % socket_matrix = randi([0 1], 6, 9)
273.             % check_ranges = [0, 1, 3, 6]
274.             % var_ranges = [0, 2, 4, 6, 9]
275.             % H = zeros(3, 4);
276.
277.             for i = 1:length(check_ranges) - 1
278.                 for j = 1:length(var_ranges) - 1
279.                     i;
280.                     j;
281.                     rows_range = check_ranges(i)+1 : + check_ranges(i+1);
282.                     cols_range = var_ranges(j)+1 : var_ranges(j+1);
283.                     submatrix = socket_matrix(rows_range, cols_range);
284.                     SUM = sum(sum(submatrix));
285.                     H(i, j) = mod(SUM, 2);
286.                 end
287.             end
288.
289.
290.             anakyklwseis = num_variable_edges - sum(sum(H));
291.             % pretty_display("Anakyklwseis: ", anakyklwseis / 2);
292.
293.         end
294.
295.         pretty_display("Anakyklwseis: ", anakyklwseis / 2);
296.         pretty_display("Tries needed: ", counter_tries);
297.         n = num_variables;
298.         k = num_variables - num_checks;
299.         rate_expected = k / n;
300.         pretty_display("Expected rate: ", rate_expected);
301.         display('*********************************************');
302.         display('*********************************************');
303.         display(' ')
304.
305.     end
306.
307.
308.
309. end
310.
311.