From 61ae4e903a9641d8027085359ef0b2229c101cec Mon Sep 17 00:00:00 2001
From: Vitaliy Filippov <vitalif@yourcmc.ru>
Date: Thu, 14 Apr 2022 20:16:09 +0300
Subject: [PATCH] WIP/experimental LRC matrix generator

---
 src/lrc/Makefile |   2 +
 src/lrc/mat.c    | 277 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 279 insertions(+)
 create mode 100644 src/lrc/Makefile
 create mode 100644 src/lrc/mat.c

diff --git a/src/lrc/Makefile b/src/lrc/Makefile
new file mode 100644
index 00000000..7c09a79f
--- /dev/null
+++ b/src/lrc/Makefile
@@ -0,0 +1,2 @@
+mat: mat.c
+	gcc -O3 -I/usr/include/jerasure -o mat mat.c -lJerasure
diff --git a/src/lrc/mat.c b/src/lrc/mat.c
new file mode 100644
index 00000000..acd62657
--- /dev/null
+++ b/src/lrc/mat.c
@@ -0,0 +1,277 @@
+#include <jerasure/reed_sol.h>
+#include <jerasure.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+// Generate LRC matrix: (groups*local + global) code rows with (data_drives) columns
+// w should be >= log2(data_drives + groups*local + global), but not necessary 8/16/32
+int* reed_sol_vandermonde_lrc_matrix(int data_drives, int groups, int local, int global, int w)
+{
+    if (w < 0 || w > 32 || data_drives + groups*local + global > (1<<w))
+    {
+        return NULL;
+    }
+    int *lrc_matrix = (int*)malloc(sizeof(int) * (local*groups+global));
+    int *matrix = reed_sol_vandermonde_coding_matrix(data_drives, local+global, w);
+    for (int gr = 0; gr < groups; gr++)
+    {
+        for (int l = 0; l < local; l++)
+        {
+            for (int j = 0; j < data_drives; j++)
+            {
+                lrc_matrix[(gr*local+l)*data_drives + j] = (j / (data_drives/groups)) == gr ? matrix[l*data_drives + j] : 0;
+            }
+        }
+    }
+    for (int i = 0; i < global; i++)
+    {
+        for (int j = 0; j < data_drives; j++)
+        {
+            lrc_matrix[(groups*local+i)*data_drives + j] = matrix[(local+i)*data_drives + j];
+        }
+    }
+    free(matrix);
+    return lrc_matrix;
+}
+
+// Check if the generated LRC with given parameters is Maximally Reconstructible (MR-LRC)
+// Example of a MR-LRC: (8, 2, 1, 2, 6, 8)
+void check_mr_lrc(int data_drives, int groups, int local, int global, int matrix_w, int w, int print)
+{
+    
+}
+
+int main()
+{
+    int W = 8, MATRIX_W = 6;
+    int n = 8, groups = 2, local = 1, global = 2;
+    //n = 4, groups = 2, local = 1, global = 1;
+    int total_rows = n+groups*local+global;
+    int *matrix = reed_sol_vandermonde_lrc_matrix(n, groups, local, global, MATRIX_W);
+    int *lrc_matrix = (int*)malloc(sizeof(int) * total_rows*n);
+    // Fill identity+LRC matrix
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < n; j++)
+            lrc_matrix[i*n + j] = j == i ? 1 : 0;
+    memcpy(lrc_matrix + n*n, matrix, (total_rows-n)*n*sizeof(int));
+    free(matrix);
+    matrix = NULL;
+    // Print LRC matrix
+    for (int i = 0; i < total_rows; i++)
+    {
+        for (int j = 0; j < n; j++)
+        {
+            printf("%d ", lrc_matrix[i*n+j]);
+        }
+        printf("\n");
+    }
+    int impossible = 0, success = 0, failures = 0;
+    int *lost_per_group = (int*)malloc(sizeof(int) * groups);
+    for (int lost = local+global+1; lost <= groups*local+global; lost++)
+    //int lost = groups*local+global;
+    {
+        int *erased_matrix = (int*)malloc(sizeof(int) * (total_rows-lost)*n);
+        int *inverted_matrix = (int*)malloc(sizeof(int) * (total_rows-lost)*n);
+        int *p = (int*)malloc(sizeof(int) * (total_rows-lost));
+        for (int i = 0; i < n; i++)
+            p[i] = i;
+        int *p2 = (int*)malloc(sizeof(int) * n);
+        if (total_rows-lost > n)
+        {
+            p[n-1] = n; // skip combinations with all N data disks (0..n-1)
+            for (int i = n; i < total_rows-lost; i++)
+                p[i] = i+1;
+            p[total_rows-lost-1]--; // will be incremented on the first step
+        }
+        int inc = total_rows-lost-1;
+        while (1)
+        {
+            p[inc]++;
+            if (p[inc] >= n+groups*local+global)
+            {
+                if (inc == 0)
+                    break;
+                inc--;
+            }
+            else if (inc+1 < total_rows-lost)
+            {
+                p[inc+1] = p[inc];
+                inc++;
+            }
+            else
+            {
+                // Check if it should be recoverable
+                for (int gr = 0; gr < groups; gr++)
+                {
+                    lost_per_group[gr] = ((gr+1)*(n/groups) > n ? (n - gr*(n/groups)) : n/groups);
+                }
+                // Calculate count of data chunks lost in each group
+                for (int j = 0; j < total_rows-lost; j++)
+                {
+                    if (j < n)
+                    {
+                        lost_per_group[(p[j] / (n/groups))]--;
+                    }
+                }
+                // Every local parity chunk is supposed to restore 1 missing chunk inside its group
+                // So, subtract local parity chunk counts from each group lost chunk count
+                for (int j = 0; j < total_rows-lost; j++)
+                {
+                    if (p[j] >= n && p[j] < n+groups*local && lost_per_group[(p[j]-n)/local] > 0)
+                    {
+                        lost_per_group[(p[j]-n)/local]--;
+                    }
+                }
+                // Every global parity chunk is supposed to restore 1 chunk of all that are still missing
+                int still_missing = 0;
+                for (int gr = 0; gr < groups; gr++)
+                {
+                    still_missing += lost_per_group[gr];
+                }
+                for (int j = 0; j < total_rows-lost; j++)
+                {
+                    if (p[j] >= n+groups*local && still_missing > 0)
+                    {
+                        still_missing--;
+                    }
+                }
+                if (still_missing <= 0)
+                {
+                    // We hope it can be recoverable. Try to invert it
+                    int invert_ok = -1;
+                    if (total_rows-lost == n)
+                    {
+                        for (int i = 0; i < n; i++)
+                            for (int j = 0; j < n; j++)
+                                erased_matrix[i*n+j] = lrc_matrix[p[i]*n+j];
+                        invert_ok = jerasure_invert_matrix(erased_matrix, inverted_matrix, n, W);
+                    }
+                    else
+                    {
+                        // Check submatrices
+                        for (int i = 0; i < n; i++)
+                            p2[i] = i;
+                        p2[n-1]--;
+                        int inc2 = n-1;
+                        while (1)
+                        {
+                            p2[inc2]++;
+                            if (p2[inc2] >= total_rows-lost)
+                            {
+                                if (inc2 == 0)
+                                    break;
+                                inc2--;
+                            }
+                            else if (inc2+1 < n)
+                            {
+                                p2[inc2+1] = p2[inc2];
+                                inc2++;
+                            }
+                            else
+                            {
+                                for (int i = 0; i < n; i++)
+                                    for (int j = 0; j < n; j++)
+                                        erased_matrix[i*n+j] = lrc_matrix[p[p2[i]]*n+j];
+                                invert_ok = jerasure_invert_matrix(erased_matrix, inverted_matrix, n, W);
+                                if (invert_ok == 0)
+                                    break;
+                            }
+                        }
+                    }
+                    if (invert_ok < 0)
+                    {
+                        failures++;
+                        printf("\nFAIL: ");
+                        for (int i = 0; i < total_rows-lost; i++)
+                        {
+                            printf("%d ", p[i]);
+                        }
+                        printf("\nDIRECT:\n");
+                        for (int i = 0; i < total_rows-lost; i++)
+                        {
+                            for (int j = 0; j < n; j++)
+                                printf("%d ", lrc_matrix[p[i]*n+j]);
+                            printf("\n");
+                        }
+                        printf("INVERSE:\n");
+                        for (int i = 0; i < total_rows-lost; i++)
+                        {
+                            for (int j = 0; j < n; j++)
+                                printf("%d ", inverted_matrix[i*n+j]);
+                            printf("\n");
+                        }
+                    }
+                    else
+                    {
+                        success++;
+                        printf("OK: ");
+                        for (int i = 0; i < total_rows-lost; i++)
+                        {
+                            printf("%d ", p[i]);
+                        }
+                        printf("\n");
+                    }
+                }
+                else
+                {
+                    impossible++;
+                    printf("IMPOSSIBLE: ");
+                    for (int i = 0; i < total_rows-lost; i++)
+                    {
+                        printf("%d ", p[i]);
+                    }
+                    printf("\n");
+                }
+            }
+        }
+        free(p2);
+        free(p);
+        free(inverted_matrix);
+        free(erased_matrix);
+    }
+    free(lost_per_group);
+    printf("\n%d recovered, %d impossible, %d failures\n", success, impossible, failures);
+    return 0;
+}
+
+// 1 1 1 1 0 0 0 0
+// 0 0 0 0 1 1 1 1
+// 1 55 39 73 84 181 225 217
+// 1 172 70 235 143 34 200 101
+//
+// Can't recover
+// 1 2 4 5 8 9 10 11    -1
+// 2 3 4 6 8 9 10 11    -1
+// FULL:
+// 1 0 0 0 0 0 0 0 
+// 0 1 0 0 0 0 0 0 
+// 0 0 1 0 0 0 0 0 
+// 0 0 0 1 0 0 0 0 
+// 0 0 0 0 1 0 0 0 
+// 0 0 0 0 0 1 0 0 
+// 0 0 0 0 0 0 1 0 
+// 0 0 0 0 0 0 0 1 
+// 1 1 1 1 0 0 0 0 
+// 0 0 0 0 1 1 1 1 
+// 1 55 39 73 84 181 225 217 
+// 1 172 70 235 143 34 200 101 
+// FIRST UNRECOVERABLE:
+// 0 1 0 0 0 0 0 0 
+// 0 0 1 0 0 0 0 0 
+// 0 0 0 0 1 0 0 0 
+// 0 0 0 0 0 1 0 0 
+// 1 1 1 1 0 0 0 0 
+// 0 0 0 0 1 1 1 1 
+// 1 55 39 73 84 181 225 217 
+// 1 172 70 235 143 34 200 101 
+// SECOND UNRECOVERABLE:
+// 0 0 1 0 0 0 0 0 
+// 0 0 0 1 0 0 0 0 
+// 0 0 0 0 1 0 0 0 
+// 0 0 0 0 0 0 1 0 
+// 1 1 1 1 0 0 0 0 
+// 0 0 0 0 1 1 1 1 
+// 1 55 39 73 84 181 225 217 
+// 1 172 70 235 143 34 200 101 
+// Ho ho ho