Decomposition

clc
clear all

%% Estw oti brhkame ta PANTA SWSTA
l_max = 3;
r_max = 6;
l_optimum = [0.3, 0.5, 0.2];
r_optimum = [0.2, 0.2 0.1, 0.1, 0.3, 0.1];
% Codewords length
n = 200;


%% Create Pi dia i and Lamda i dia i
r_i_dia_i_list = zeros(1, r_max);
l_i_dia_i_list = zeros(1, l_max);
for index = 1:r_max
    r_i_dia_i_list(index) = r_optimum(index) / (index + 1);     % Epeidi arxizei apo to 2 enw emeis apo to 1
end

for index = 1:l_max
    l_i_dia_i_list(index) = l_optimum(index) / (index + 1);
end
SUM = sum(l_i_dia_i_list);

r_i_dia_i_list;
l_i_dia_i_list;

%% Create LAMDA AND RHO LISTS
RHO = r_i_dia_i_list / SUM * n;
LAMDA = l_i_dia_i_list / SUM * n;
RHO_HAT = floor(RHO);
LAMDA_HAT = floor(LAMDA);
% Calculate A
A_condition = sum(RHO - RHO_HAT);


%% GENERAL CONSTRAINTS
% Declare 2 sum constant variables
constant_19 = n - sum(LAMDA_HAT);
i_r = 2:r_max+1;
i_l = 2:l_max+1;
constant_20 = sum(i_r .* RHO_HAT) - sum(i_l .* LAMDA_HAT);

% Prwta stoixeia ta "l", kai meta ta "r"
Aeq = zeros(2, l_max + r_max);
Beq = zeros(2, 1);

% Constraint 19
Aeq(1, 1:l_max) = 1;
Beq(1, 1) = constant_19;
% Constraint 20
Aeq(2, 1:l_max) = i_l;
Aeq(2, l_max+1 : l_max+r_max) = - i_r;
Beq(2, 1) = constant_20;
% Constraint 21


%% CASE CONSTRAINTS
% CASE 1 -  ANISOTIKOS - MATRIX A1
A1 = zeros(1, l_max + r_max);
A1(1, l_max+1 : l_max+r_max) = -1;
B1 = -ceil(A_condition);

% CASE 2
A2 = zeros(1, l_max + r_max);
A2(1, l_max+1 : l_max+r_max) = 1;
B2 = floor(A_condition);

% OBJECTIVE FUNCTIONS
% CASE 1
f1 = zeros(1, l_max + r_max);
f1(1, l_max+1 : l_max+r_max) = 1;

% CASE 2
f2 = -f1;


%% SIMULATIONS
lb = zeros(1, l_max + r_max);
ub = ones(1, l_max + r_max);
% CASE 1
result1 = linprog(f1, A1, B1 , Aeq, Beq, lb, ub);
% CASE 2
result2 = linprog(f2, A2, B2 , Aeq, Beq, lb, ub);

%% FIND THE BEST
error1 = (sum(result1(l_max+1 : l_max+r_max)) - A_condition) ^ 2;
error2 = (sum(result2(l_max+1 : l_max+r_max)) - A_condition) ^ 2;
if error1 >= error2
    result = result2;
else
    result = result1;
end
% result = [xl1, xl2, xl3, xr1, ..., xr6] ---> length = 9


%% FIND THE NEW LAMDA AND RHO - UPDATE
n;
result = result';
LAMDA_HAT;
RHO_HAT;
LAMDA_LIST = LAMDA_HAT + result(1 : l_max);
RHO_LIST = RHO_HAT + result(l_max + 1 : l_max + r_max);


%% DISPLAY AKMES
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% MPOREI KAPOIA NOUMERA NA MIN EINAI INTEGERS, ENW ANAPARISTOUN # AKMWN !
num_variable_edges = sum(LAMDA_LIST .* i_l);
num_check_edges = sum(RHO_LIST .* i_r);
num_variables = sum(LAMDA_LIST);
num_checks = sum(RHO_LIST);
%   Above: L'(1) = P'(1)  (sxesh 2 kai 3)

pretty_display("n = ", n);
pretty_display("Variable degrees = ", LAMDA_LIST);
pretty_display("Check degrees = ", RHO_LIST);
pretty_display("Num of variable edges = ", num_variable_edges);
pretty_display("Num of check edges = ", num_check_edges);
pretty_display("Num of VARIABLES = ", num_variables);
pretty_display("Num of CHECKS = ", num_checks);



%% MATRIX 545 x 545 (about sockets)
sockets = num_check_edges;
checks = num_checks;
vars = num_variables;

% Socket Matrix is 545 x 545
socket_matrix = zeros(sockets);
comb1 = 1:sockets;
comb2 = comb1(randperm(length(comb1)));

for i = 1:sockets
    ch = comb1(i);
    var = comb2(i);
    socket_matrix(ch, var) = 1;
end


%% Create a vector with 82 2's and 91 3's and 27 4's (n=200)
var_temp = zeros(1, vars);
counter = 1;
for i = 1:length(LAMDA_LIST)
    var_temp(counter : counter + LAMDA_LIST(i) - 1) = ...
    (i+1) * ones(1, LAMDA_LIST(i));
    counter = LAMDA_LIST(i) + counter;
end
var_temp;


check_temp = zeros(1, checks);
counter = 1;
for i = 1:length(RHO_LIST)
    check_temp(counter : counter + RHO_LIST(i) - 1) = ...
    (i+1) * ones(1, RHO_LIST(i));
    counter = RHO_LIST(i) + counter;
end
check_temp;


%% H is a decomposition of socket matrix
H = zeros(checks, vars);
% for i = 1:vars
%     for j = 1:checks
%         MAT = sockets( : , :);
%         H(i, j) = sum(MAT);
%     end
% end




%% ACCUMULATIVE RANGES
var_ranges = zeros(1, length(var_temp));
var_ranges = [0 var_ranges];

for i = 1:length(var_temp)
    S = var_ranges(i) + var_temp(i);
    var_ranges(i+1) = S;
end
var_ranges;


check_ranges = zeros(1, length(check_temp));
check_ranges = [0 check_ranges];

for i = 1:length(check_temp)
    S = check_ranges(i) + check_temp(i);
    check_ranges(i+1) = S;
end
check_ranges;


%% Decomposition onto submatrices
% var_ranges and check-ranges have 1 more element (0 in the beginning)
socket_matrix = randi([0 1], 6, 9)
check_ranges = [0, 1, 3, 6]
var_ranges = [0, 2, 4, 6, 9]
H = zeros(3, 4);

for i = 1:length(check_ranges) - 1
    for j = 1:length(var_ranges) - 1
        i
        j
        rows_range = check_ranges(i)+1 : + check_ranges(i+1)
        cols_range = var_ranges(j)+1 : + var_ranges(j+1)
        submatrix = socket_matrix(rows_range, cols_range)
        SUM = sum(sum(submatrix))
        H(i, j) = mod(SUM, 2)
    end
end

H