Advertisement
makispaiktis

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.  
Advertisement
Advertisement
Advertisement
RAW Paste Data Copied
Advertisement