Avalanche RLP encoded block translated into matrix format and projected into SL(667*2, 16)

avax_rlp_matrix_form.cpp

#include <iostream>
#include <map>
#include <vector>
#include <string> 
#include <cstdint>

#include <armadillo>

#include <sys/stat.h>
#include <json/json.h>
#include <json/reader.h>

std::string alpha = "0123456789abcdef";

/*
    
    clang++ -fcolor-diagnostics block_encoding_text.cpp -o bt -Wall -O3 -std=gnu++14 -stdlib=libc++ -pedantic -I/usr/local/include `pkg-config --cflags --libs jsoncpp` /usr/local/lib/libarmadillo.dylib

    curl -X POST --data '{ "jsonrpc":"2.0", "id" :1, "method" :"debug_getBlockRlp", "params" : [block_num]}' -H 'content-type:application/json;' http://127.0.0.1:9545/ext/bc/C/rpc | python -m json.tool > block_num.json ; ./b1 "block_num.json"

    time ./bt "1.json" "10.json" "1.0"
*/

const int buff_size = 8;

// Lookup table of Fibonacci numbers that can fit in a ulong.
const static std::vector<uint64_t> fibLookup{
    1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657,
    46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352,
    24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903,
    2971215073, 4807526976, 7778742049, 12586269025, 20365011074, 32951280099, 53316291173, 86267571272, 
    139583862445, 225851433717, 365435296162, 591286729879, 956722026041, 1548008755920, 2504730781961, 
    4052739537881, 6557470319842, 10610209857723, 17167680177565, 27777890035288, 44945570212853, 
    72723460248141, 117669030460994, 190392490709135, 308061521170129, 498454011879264, 806515533049393, 
    1304969544928657, 2111485077978050, 3416454622906707, 5527939700884757, 8944394323791464, 14472334024676221,
    23416728348467685, 37889062373143906, 61305790721611591, 99194853094755497, 160500643816367088, 
    259695496911122585, 420196140727489673, 679891637638612258, 1100087778366101931, 1779979416004714189, 
    2880067194370816120, 4660046610375530309, 7540113804746346429, 12200160415121876738ULL
};

//https://stackoverflow.com/questions/3543783/what-is-the-c-analog-of-c-sharp-byte c# byte -> c++ uint8_t

/// The most significant bit in a byte.
const uint8_t fibMSB = 0x80; //8 bytes -> 128 or 0x80, 16 bytes -> 32768 or 0x8000

/// Maximum possible length of an encoded symbol.
const int MAX_ENCODED_LENGTH = buff_size;

void fib_encode_many(std::vector<unsigned short>& values, int offset=0, int count=0) {
    // Allocate working buffer the max length of an encoded UInt64 + 1 byte
    std::ofstream stream("instructions.bin.fib", std::ios::binary);
    uint8_t buffer[MAX_ENCODED_LENGTH];
    int bitOffset = 0;

    // need to initialize inorder to avoid weird bit flipping behavoir
    for (int i=0; i< MAX_ENCODED_LENGTH; i++){
        buffer[i] = (uint8_t)0;
    }

    // Iterate through all symbols
    for (int valueIdx = offset; valueIdx < offset + count; valueIdx++) {
        unsigned short value = values[valueIdx];
        int residualBits = 0;
        // Reset size for next symbol
        int maxByte = -1;
        
        // Fibonacci doesn't support 0s, so add 1 to allow for them
        value++;

        // #1 Find the largest Fibonacci number equal to or less than N; subtract this number from N, keeping track of the remainder.
        // #3 Repeat the previous steps, substituting the remainder for N, until a remainder of 0 is reached.
        for (int fibIdx = fibLookup.size() - 1; fibIdx >= 0; fibIdx--) {
            // #2 If the number subtracted was the ith Fibonacci number F(i), put a 1 in place i−2 in the code word (counting the left most digit as place 0).
            if (value >= fibLookup[fibIdx]) {
                // Calculate offsets
                int adjustedIdx = fibIdx + bitOffset;
                int byteIdx = adjustedIdx / buff_size;
                int bitIdx = adjustedIdx % buff_size;

                // If this is the termination fib, add termination bit
                if (-1 == maxByte) {
                    // Note parameters for this symbol
                    maxByte = (adjustedIdx + 1) / buff_size;
                    residualBits = (adjustedIdx + 2) % buff_size; // Add two bits being written

                    // Append bits to output
                    int terminationByteIdx = (adjustedIdx + 1) / buff_size;
                    int terminationBitIdx = (adjustedIdx + 1) % buff_size;
                    buffer[terminationByteIdx] |= (uint8_t)(0x01 << (buff_size - 1 - terminationBitIdx)); // Termination bit
                }
                // Flag that fib is used
                buffer[byteIdx] |= (uint8_t)(0x01 << (buff_size - 1 - bitIdx)); // Flag bit
                // Deduct Fibonacci number from value
                value -= fibLookup[fibIdx];
            }
        }

        // Write n-1 output bytes
        for (int j = 0; j < maxByte; j++) {
            stream.write(reinterpret_cast<const char *>(&buffer[j]), sizeof(uint8_t));
            buffer[j] = 0;
        }

        // Write last byte if complete, or no more symbols to encode
        if (residualBits == 0 || // No residual bits
            valueIdx == offset + count - 1 // No more symbols to process
        ) {
            stream.write(reinterpret_cast<const char *>(&buffer[maxByte]), sizeof(uint8_t));
            buffer[maxByte] = 0;
        } else if (maxByte > 0) {
            buffer[0] = buffer[maxByte];
            buffer[maxByte] = 0;
        }

        bitOffset = residualBits;
    }

    std::cout << "\t\t" << stream.tellp() << std::endl;
    stream.close();
}

std::vector<unsigned short> fib_decode_many(int offset=0) {
    // Current symbol being decoded
    std::vector<unsigned short> values;
    uint64_t symbol = 0;

    // Next Fibonacci number to test
    int nextFibIndex = 0;

    // State of the last bit while decoding
    bool lastBit = false;

    /*
        TODO: read to stringstream and decompress
    */
    std::ifstream ifd;
    ifd.open("instructions.bin.fib", std::ios::binary | std::ios::ate);
    int size = ifd.tellg();
    int total_intructions = size / sizeof (uint8_t);
    ifd.seekg(0,  std::ios::beg);
    for (int i=0; i < total_intructions; i++) {
        // Read byte of input, and throw error if unavailable
        uint8_t b;
        ifd.read(reinterpret_cast<char *>(&b), sizeof(uint8_t));


        if (b < 0) {
            std::cout << "Input ends with a partial symbol. More bytes required to decode." << std::endl;
            exit(0);
        }

        // For each bit of buffer
        for (int bi = 0; bi < buff_size; bi++) {
            // If bit is set... MSB == 0x8000
            //uint16_t b_shifted = b << 1;
            if (((b << bi) & fibMSB) > 0) {
                // If double 1 bits
                if (lastBit) {
                    // Remove zero offset
                    symbol--;

                    if (symbol != 0) {
                        // Reset for next symbol
                        values.push_back((unsigned short) symbol);
                        symbol = 0;
                        nextFibIndex = 0;
                    } else {
                        symbol++;
                    }

                    lastBit = false;
                    continue;
                }
                // Add value to current symbol
                //unsigned short pre = symbol;
                symbol += fibLookup[nextFibIndex];
                // Note bit for next cycle
                lastBit = true;
            } else {
                // Note bit for next cycle
                lastBit = false;
            }

            // Increment bit position
            nextFibIndex++;
        }
    }
    ifd.close();
    return values;
}

Json::Value load_json(std::string& file){
    Json::Value root;
    {
        try{
            std::ifstream json_data(file);
            json_data >> root;
        }catch(...){
             Json::Value root;
        }
    }
    return root;
}

std::string to_hex(unsigned short& d) {
    unsigned short r = d % 16;
    if ((d - r) == 0) {
        return alpha.substr(r,1);
    }
    unsigned short nd = (d - r) / 16;
    return to_hex(nd) + alpha.substr(r,1);
}

std::string join_string(std::vector<unsigned short>& v, std::string delimiter, bool print_index=false){
    int vector_len = v.size();
    std::string ov;
    for(int i = 0; i< vector_len; i++){
        if(i < vector_len-1){
            if (print_index)
                ov += std::to_string(i) + ": " + std::to_string(v[i]) + delimiter;
            else
                ov += std::to_string(v[i]) + delimiter;
        }else{
            if (print_index)
                ov += std::to_string(i) + ": " + std::to_string(v[i]);
            else
                ov += std::to_string(v[i]);
        }
    }
    return ov;  
}

double compute_abs_det(arma::imat& input){
    //https://stackoverflow.com/questions/47135703/determinant-of-sparse-matrix
    /* avoid deterniant calc for large matrix */
    arma::dmat Q;
    arma::dmat R;
    arma::dmat d_input = arma::conv_to<arma::dmat>::from(input);

    arma::qr(Q, R, d_input);

    arma::dvec diag_R = R.diag();
    double sum_det = 1;
    for (int i =0; i< diag_R.size(); i++){
        sum_det *= std::abs(diag_R(i));
    }
    return sum_det;
}

arma::sp_imat elementary_row_add(arma::sp_imat E, arma::sp_imat& input, int i, int j, double value){
    if (i == j || input.n_rows != input.n_cols)
        return input;
    else {
        E(i,j) = value;
        return E * input;
    }
}

std::vector<arma::imat> load_block_proj_mod(){

    std::string block_data_str = alpha;
    int max_block_data_row = 667; //10000
    arma::imat block_rep(max_block_data_row, max_block_data_row);
    int max_alpha_size = block_data_str.size();
    for (int i=0; i<max_alpha_size; i++) {
        int mod_val = (i % max_alpha_size);
        block_rep(i, i) = mod_val;
    }

    // https://math.stackexchange.com/questions/2817168/can-a-non-invertible-matrix-be-extended-to-an-invertible-one
    /*
        block_rep -> non invertable
        block_rep -> proj_block_rep makes (special linear group: real matrices with determinant 1)
        https://en.wikipedia.org/wiki/Special_linear_group
        SL(667*2, 16) -> 667*2 * 16 (21,344) unique orderings of elements
        https://en.wikipedia.org/wiki/File:SL(2,3);_Cayley_table.svg for example

        https://en.wikipedia.org/wiki/Special_linear_group#Generators_and_relations -> https://en.wikipedia.org/wiki/Elementary_matrices
        how many elemntary matricies to converte proj_block_rep_1 to proj_block_rep_2?
            - row addition

        to do need to experement with computing the differences between two blocks and then for each different, which set of elemetary operations need to be done and in which order
            -> which linear combination of previous blocks rows will equal current block rows? 


    */
    arma::imat proj_block_rep(max_block_data_row * 2, max_block_data_row * 2);
    arma::imat proj_block_rep_inv(max_block_data_row * 2, max_block_data_row * 2);
    arma::imat I,I2;
    I.eye(max_block_data_row, max_block_data_row);
    I2.eye(max_block_data_row*2, max_block_data_row*2);
    // X.submat( first_row, first_col, last_row, last_col )
    proj_block_rep.submat(0,0, max_block_data_row-1, max_block_data_row-1) = block_rep;
    proj_block_rep.submat(max_block_data_row,0, (max_block_data_row*2)-1, (max_block_data_row-1)) = I;
    proj_block_rep.submat(0,max_block_data_row, max_block_data_row-1, (max_block_data_row*2)-1) = I;


    proj_block_rep_inv.submat(max_block_data_row,0, (max_block_data_row*2)-1, (max_block_data_row-1)) = I;
    proj_block_rep_inv.submat(0,max_block_data_row, max_block_data_row-1, (max_block_data_row*2)-1) = I;
    proj_block_rep_inv.submat(max_block_data_row,max_block_data_row, (max_block_data_row*2)-1, (max_block_data_row*2)-1) = -1 * block_rep;


    arma::imat c_proj_block_rep_t = arma::conj(proj_block_rep).t();
    arma::imat test_I2 = c_proj_block_rep_t * proj_block_rep;

    bool is_identity = arma::approx_equal(test_I2, I2, "absdiff", 0);

    std::cout << "block_rlp_len(artificial): " << block_data_str.size() << std::endl;
    std::cout << "is_identity(artificial): " << is_identity << std::endl;
    std::cout << "abs_det(proj_block_rep->artificial): " << compute_abs_det(proj_block_rep) << std::endl;
    std::cout << "abs_det(c_proj_block_rep_t->artificial): " << compute_abs_det(c_proj_block_rep_t) << std::endl;
    return {proj_block_rep, proj_block_rep_inv};
}

std::vector<arma::imat> load_block_proj(std::string& filename, std::map<std::string, int>& char_num_map){
    Json::Value block_data;
    struct stat buf;
    if((stat(filename.c_str(), &buf) != 0)){
        std::cout << "FAIL" << std::endl;
        exit(0);
    }

    try{
        block_data = load_json(filename);
    }catch(...){
        std::cout << "FAIL" << std::endl;
        exit(0);
    }

    std::string block_data_str = block_data["result"].asString();
    int max_block_data_row = 667; //10000
    arma::imat block_rep(max_block_data_row, max_block_data_row);
    int row_start = 0;
    for (int i=0; i<block_data_str.size(); i++) {
        if (((i % max_block_data_row) == 0) && (i > 0))
            row_start++;
        block_rep(row_start,i % max_block_data_row) = char_num_map[block_data_str.substr(i,1)];
    }

    // https://math.stackexchange.com/questions/2817168/can-a-non-invertible-matrix-be-extended-to-an-invertible-one
    /*
        block_rep -> non invertable
        block_rep -> proj_block_rep makes (special linear group: real matrices with determinant 1)
        https://en.wikipedia.org/wiki/Special_linear_group
        SL(667*2, 16) -> 667*2 * 16 (21,344) unique orderings of elements
        https://en.wikipedia.org/wiki/File:SL(2,3);_Cayley_table.svg for example

        https://en.wikipedia.org/wiki/Special_linear_group#Generators_and_relations -> https://en.wikipedia.org/wiki/Elementary_matrices
        how many elemntary matricies to converte proj_block_rep_1 to proj_block_rep_2?
            - row addition

        to do need to experement with computing the differences between two blocks and then for each different, which set of elemetary operations need to be done and in which order
            -> which linear combination of previous blocks rows will equal current block rows? 


    */
    arma::imat proj_block_rep(max_block_data_row * 2, max_block_data_row * 2);
    arma::imat proj_block_rep_inv(max_block_data_row * 2, max_block_data_row * 2);
    arma::imat I,I2;
    I.eye(max_block_data_row, max_block_data_row);
    I2.eye(max_block_data_row*2, max_block_data_row*2);
    // X.submat( first_row, first_col, last_row, last_col )
    proj_block_rep.submat(0,0, max_block_data_row-1, max_block_data_row-1) = block_rep;
    proj_block_rep.submat(max_block_data_row,0, (max_block_data_row*2)-1, (max_block_data_row-1)) = I;
    proj_block_rep.submat(0,max_block_data_row, max_block_data_row-1, (max_block_data_row*2)-1) = I;


    proj_block_rep_inv.submat(max_block_data_row,0, (max_block_data_row*2)-1, (max_block_data_row-1)) = I;
    proj_block_rep_inv.submat(0,max_block_data_row, max_block_data_row-1, (max_block_data_row*2)-1) = I;
    proj_block_rep_inv.submat(max_block_data_row,max_block_data_row, (max_block_data_row*2)-1, (max_block_data_row*2)-1) = -1 * block_rep;


    arma::imat c_proj_block_rep_t = arma::conj(proj_block_rep).t();
    arma::imat test_I2 = c_proj_block_rep_t * proj_block_rep;

    bool is_identity = arma::approx_equal(test_I2, I2, "absdiff", 0);

    std::cout << "block_rlp_len(" << filename << "): " << block_data_str.size() << std::endl;
    std::cout << "is_identity(" << filename << "): " << is_identity << std::endl;
    std::cout << "abs_det(proj_block_rep->" << filename << "): " << compute_abs_det(proj_block_rep) << std::endl;
    std::cout << "abs_det(c_proj_block_rep_t->" << filename << "): " << compute_abs_det(c_proj_block_rep_t) << std::endl;

    /*
    arma::sp_imat proj_block_rep_sparce = arma::sp_imat(proj_block_rep);
    proj_block_rep_sparce.save("proj_block_rep_sparce_" + filename + ".mat");
    proj_block_rep.save("proj_block_rep_raw_" + filename + ".mat", arma::raw_ascii);
    */

    return {proj_block_rep, proj_block_rep_inv};
}

std::map<std::string, int> init_char_map(std::string& alphabet){
    std::map<std::string, int> char_num_map;
    for(int i=0; i< alphabet.size(); i++) {
        char_num_map[alphabet.substr(i,1)] = i;
    }
    return char_num_map;
}

void write_decode_instructions(std::vector<unsigned short>& inst) {
    //std::cout << "inst.size(): " << inst.size() << std::endl;
    /*
        TODO: write to stringstream and compress
    */
    fib_encode_many(inst, 0, inst.size());
}

std::vector<unsigned short> read_decode_instructions() {
    /*
        TODO: read to stringstream and decompress
    */
    std::vector<unsigned short> instructions = fib_decode_many(0);

    int neg_offset = 16;
    int skip_offset = 32;

    int start = instructions[0];
    int end = instructions[1];
    std::cout << "start: " << start << std::endl;
    std::cout << "end: " << end << std::endl;

    arma::sp_imat E;
    E.eye(667*2, 667*2);
    arma::sp_imat a1_sparce = E;
    int current_row = 0;
    std::map<int, std::vector<int>> row_column_intstructions;
    std::map<int, int> row_col_index;
    int total_intructions = instructions.size();
    for (int i=2; i < total_intructions; i++){
        if ((instructions[i] > 15+neg_offset) && (instructions[i] < skip_offset+neg_offset)) {
            if (instructions[i+1] != neg_offset){
                //std::cout << "count in row: " << instructions[i] << std::endl;
                if (current_row == 0)
                    current_row = instructions[0];
                else
                    current_row++;

                row_col_index[current_row] = 0;
            } else {
                //std::cout << "\t\tskip columns: "<< instructions[i] << std::endl;
                if (instructions[i] != neg_offset) {
                    row_column_intstructions[current_row].push_back(instructions[i]);

                    if ((instructions[i] < 16+neg_offset)){
                        // need to take the opposite to undo transformation
                        int op_val = (instructions[i] - neg_offset) * -1;
                        int x = current_row;
                        int c = row_col_index[x];
                        a1_sparce = elementary_row_add(E, a1_sparce, x, (c == x) ? (x - 1) : c, op_val);
                        row_col_index[current_row]++;
                    } else {
                        // skip value
                        int skip_value = instructions[i] - neg_offset - skip_offset;
                        row_col_index[current_row] += skip_value;
                    }
                }
            }
        } else {
            row_column_intstructions[current_row].push_back(instructions[i]);

            if ((instructions[i] < 16+neg_offset)){
                // need to take the opposite to undo transformation
                int op_val = (instructions[i] - neg_offset) * -1;
                int x = current_row;
                int c = row_col_index[x];
                a1_sparce = elementary_row_add(E, a1_sparce, x, (c == x) ? (x - 1) : c, op_val);
                row_col_index[current_row]++;
            } else {
                // skip value
                int skip_value = instructions[i] - neg_offset - skip_offset;
                row_col_index[current_row] += skip_value;
            }
        }
    }

    arma::imat a = arma::imat(a1_sparce);
    a.save("a_b0->b1_undo.mat", arma::raw_ascii);
    return instructions;
}

int main(int argc, char *argv[]){

    std::string sample_block = (std::string)argv[1];
    std::string sample_block_1 = (std::string)argv[2];

    std::cout << "sample_block: " << sample_block << std::endl;
    std::cout << "sample_block_1: " << sample_block_1 << std::endl;
    int neg_offset = 16;
    int skip_offset = 32;


    if (sample_block != "" && sample_block_1 != "") {

        std::map<std::string, int> char_num_map = init_char_map(alpha);
        std::vector<arma::imat> b0  = load_block_proj_mod();
        std::vector<arma::imat> b1 = load_block_proj(sample_block_1, char_num_map);

        /* 
            B1 = B0 * A 
            B0^-1 * B1 = A
        */

        arma::imat a = b0[1] * b1[0];
        a.save("a_b0->b1.mat", arma::raw_ascii);
        int non_one_sum_count = 0;
        int emra_iter_count = 0;
        int emra_mod_count = 0;
        int prev_blank_count = 0;


        arma::sp_imat E;
        E.eye(a.n_rows, a.n_cols);
        arma::sp_imat a_sparce = arma::sp_imat(a);
        arma::sp_imat a1_sparce = a_sparce;
        int start_offset = 0;
        int last_offset = 0;
        std::string decoding_string = "";
        std::vector<unsigned short> all_int_operations;
        int row_mod_idx = 0;
        std::map<int, std::vector<unsigned short>> row_ops_map;
        for (int i=0; i< a.n_rows; i++) {
            double sum = arma::accu(a.row(i));
            if (sum != 1.0) {
                if (start_offset == 0)
                    start_offset = i;
                
                non_one_sum_count++;
                last_offset = i;
                int prev_blank = 0;

                /*
                arma::uvec a_nz = arma::find(a.row(i) > 0);
                std::cout << "a.row("<< i <<", " << a_nz.size()<<") sum:" << sum << std::endl;
                */
                int row_mod_count = 0;
                std::string operations = "";
                std::vector<unsigned short> int_operations;
                std::vector<unsigned short> int_operations_pure;

                
                for (int j = 0; j <= i;  j++) {
                    /*
                        need to determine how many elementary_row_add need to happen

                        for each row, need to keep track of
                            - next non blank col
                            - total operations per row
                            - operation value

                        TODO: how to minimize the number of operations

                    */

                    int val = a1_sparce(i, j);
                    if (val != 0) {
                        if (prev_blank != 0){
                            // these set of instruction tells you that you should skip x many cells before starting elementary operations
                            //unsigned short tmp = (unsigned short)(prev_blank + neg_offset + skip_offset);
                            int_operations.push_back(prev_blank + neg_offset + skip_offset);
                            //int_operations_pure.push_back((unsigned short)(prev_blank + neg_offset));
                            prev_blank_count++;
                        }
                        int op_val;
                        int mod_row;
                        if (j < i) {
                            op_val = val * -1;
                            mod_row = j;
                        } else {
                            op_val = (val * -1) + 1;
                            mod_row = i-1;
                        }
                        a1_sparce = elementary_row_add(E, a1_sparce, i, mod_row, op_val);
                        //unsigned short tmp1 = (unsigned short)(op_val+neg_offset);
                        int_operations.push_back(op_val+neg_offset);
                        //int_operations_pure.push_back(tmp1);
                        emra_mod_count++;
                        row_mod_count++;
                        prev_blank = 0;
                    } else {
                        prev_blank++;
                    }
                    emra_iter_count++;
                    
                }
                row_mod_idx++;
                all_int_operations.push_back(16+neg_offset);
                all_int_operations.insert(all_int_operations.end(), int_operations.begin(), int_operations.end());
            }
            
        }

        std::vector<unsigned short> o_vec;
        o_vec.resize(2+all_int_operations.size());
        o_vec.push_back(((unsigned short)start_offset));
        o_vec.push_back(((unsigned short)last_offset));
        o_vec.insert(o_vec.end(), all_int_operations.begin(), all_int_operations.end());

        write_decode_instructions(o_vec);

        std::cout << "non_one_sum_count:" << non_one_sum_count << std::endl;
        std::cout << "emra_iter_count:" << emra_iter_count << std::endl;
        std::cout << "emra_mod_count:" << emra_mod_count << std::endl;
        std::cout << "prev_blank_count:" << prev_blank_count << std::endl;
        std::cout << "start_offset:" << start_offset << std::endl;
        std::cout << "last_offset:" << last_offset << std::endl;
        std::cout << "total operations:" << emra_mod_count + prev_blank_count + (last_offset - start_offset) << std::endl;
        
        a_sparce.save("a_b0->b1_sparce.mat");
        arma::imat a1 = arma::imat(a1_sparce);
        a1.save("a1_b0->b1.mat", arma::raw_ascii);
        bool is_identity = arma::approx_equal(a1, arma::imat(E), "absdiff", 0);
        std::cout << "is_identity(a1_b0->b1.mat): " << is_identity << std::endl;
    } else {
        std::vector<unsigned short> instructions = read_decode_instructions();
    }
}