host_utils.h

Go to the documentation of this file.

#ifndef HOST_H
#define HOST_H

#include <string>
#include <vector>
#include <array>
#include <complex>
#include <random>
#include <stack>
#include <map>
#include <iostream>

using namespace std;

#define AL_END_H 0b000  // 0 stopping condition
#define AL_INS_H 0b001  // 1 Align Insertion
#define AL_MMI_H 0b010  // 2 Align Match/Mismatch
#define AL_DEL_H 0b011  // 3 Align Deletion
#define AL_NULL_H 0b100 // 4 Do not change coordinate

namespace Random
{
    // Generate a random sequence of length `length` with alphabet size `alphabet_size`.
    // The alphabet is [0, alphabet_size).
    template <int N>
    string Sequence(char alphabet[N], int length)
    {
        std::random_device rd;

        // Use the seed from random_device to seed the random engine
        std::mt19937 rng(rd());

        // Create a uniform distribution for random numbers between 1 and 100
        std::uniform_int_distribution<int> dist(1, 100);

        string seq = "";
        for (int i = 0; i < length; i++)
        {
            seq += alphabet[dist(rng) % N];
        }
        return seq;
    }

    template <typename T>
    std::vector<T> SequenceComplex(int length)
    {
        std::random_device rd;

        // Use the seed from random_device to seed the random engine
        std::mt19937 rng(rd());

        // Create a uniform distribution for random numbers between 1 and 100
        std::uniform_int_distribution<int> dist(1, 100);

        std::vector<T> seq;
        for (int i = 0; i < length; i++)
        {
            seq.push_back(T(dist(rng) % 10, dist(rng) % 10));
        }
        return seq;
    }

    template <int WID>
    std::vector<std::array<int, 5>> SequenceProfileAlignment(int len)
    {
        std::vector<std::array<int, 5>> profile;
        for (int i = 0; i < len; i++)
        {
            int remaining = WID;
            std::array<int, 5> col;
            for (int j = 0; j < 4; j++)
            {
                col[j] = rand() % remaining;
                remaining -= col[j];
            }
            col[4] = WID - remaining;
            profile.push_back(col);
        }
        return profile;
    }

}

namespace HostUtils
{
    namespace IO
    {
        map<string, std::vector<string>> read_sequences_from_json(string file_path);

        std::vector<std::string> readFasta(const std::string &filePath);

        template <typename T, size_t M, size_t N>
        void SwitchDimension(T (&src)[M][N], T (&dst)[N][M])
        {
            for (int i = 0; i < M; i++)
            {
                for (int j = 0; j < N; j++)
                {
                    dst[j][i] = src[i][j];
                }
            }
        }
    }

    namespace Generate
    {
    }

    namespace Sequence
    {
        int base_to_num(char base);

        char num_to_base(int num);

        char num_to_aa(int num);
        int aa_to_num(char aa);

        template <typename T>
        string nav_to_string(T nav)
        {
            // Copilot, please write the above switch statement in to if statement
            if (nav == (T) AL_NULL_H)
            {
                return "";
            }
            else if (nav == (T) AL_MMI_H)
            {
                return "D";
            }
            else if (nav == (T) AL_INS_H)
            {
                return "L";
            }
            else if (nav == (T) AL_DEL_H)
            {
                return "U";
            }
            else if (nav == (T) AL_END_H)
            {
                return "*";
            }
            else
            {
                return "?";
            }
        }

        template <typename T, size_t MAX_QRY_LENGTH, size_t MAX_REF_LENGTH>
        std::vector<string> ReconstructTracebackProfile(std::vector<string> querys, std::vector<string> references,
                                                        int query_start_idx, int reference_start_idx,
                                                        T (&tb_streams)[MAX_REF_LENGTH + MAX_QRY_LENGTH])
        {
            std::vector<string> alignments_query;
            std::vector<string> alignments_reference;
            // print size of query
            int num_queries = querys.size();
            int num_references = references.size();
            // initialize alignments
            for (int i = 0; i < num_queries; i++)
            {
                alignments_query.push_back("");
            }
            for (int i = 0; i < num_references; i++)
            {
                alignments_reference.push_back("");
            }
            std::cout << "Number of queries: " << num_queries << endl;
            std::cout << "Number of references: " << num_references << endl;
            // If we also want to match the portion for the alignment exceeds the start index, use this
            // string alignment_query = query.substr(query_start_idx + 1, query.length());
            // string alignment_reference = reference.substr(reference_start_idx + 1, reference.length());

            std::vector<stack<char>> query_stack;
            std::vector<stack<char>> reference_stack;

            // Initialize vectors for querys
            for (int i = 0; i < querys.size(); i++)
            {
                stack<char> q;
                for (int j = 0; j < query_start_idx; j++)
                {
                    q.push(querys[i][j]);
                }
                query_stack.push_back(q);
            }
            // Initialize vectors for references
            for (int i = 0; i < references.size(); i++)
            {
                stack<char> r;
                for (int j = 0; j < reference_start_idx; j++)
                {
                    r.push(references[i][j]);
                }
                reference_stack.push_back(r);
            }

            // print traceback stream
            // printf("Traceback stream: ");
            // for (int i = 0; i < MAX_REFERENCE_LENGTH + MAX_QUERY_LENGTH; i++) {
            //     printf("%d ", tb_stream[i].to_int());
            // }

            T *curr_ptr = &tb_streams[0];
            // Insert the characters from the stack to the alignment strings to their beginning.
            // Iterating in order the tb_steram
            int cnt = 0;
            while (*curr_ptr != (T)AL_END_H)
            {
                // printf("curr_ptr: %d\n", curr_ptr->to_int());
                if (*curr_ptr == (T)AL_MMI_H)
                {
                    for (int i = 0; i < querys.size(); i++)
                    {
                        alignments_query[i].insert(0, 1, query_stack[i].top());
                        query_stack[i].pop();
                    }
                    for (int i = 0; i < reference_stack.size(); i++)
                    {
                        alignments_reference[i].insert(0, 1, reference_stack[i].top());
                        reference_stack[i].pop();
                    }
                }
                else if (*curr_ptr == (T)AL_INS_H)
                {
                    for (int i = 0; i < query_stack.size(); i++)
                    {
                        alignments_query[i].insert(0, 1, '_');
                    }
                    for (int i = 0; i < reference_stack.size(); i++)
                    {
                        alignments_reference[i].insert(0, 1, reference_stack[i].top());
                        reference_stack[i].pop();
                    }
                }
                else if (*curr_ptr == (T)AL_DEL_H)
                {
                    for (int i = 0; i < query_stack.size(); i++)
                    {
                        alignments_query[i].insert(0, 1, query_stack[i].top());
                        query_stack[i].pop();
                    }
                    for (int i = 0; i < reference_stack.size(); i++)
                    {
                        alignments_reference[i].insert(0, 1, '_');
                    }
                }
                else if (*curr_ptr == (T)AL_NULL_H)
                {
                    // Do nothing, AL_NULL Doesn't change the position
                }
                else
                {
                    printf("Alignment Output Iteartion End\n");
                }
                curr_ptr++;
                cnt++;
            }

            // Finishing up concatenating the rest of the characters in the stack

            for (int i = 0; i < query_stack.size(); i++)
            {
                while (!query_stack[i].empty())
                {
                    alignments_query[i].insert(0, 1, query_stack[i].top());
                    query_stack[i].pop();
                }
            }

            for (int i = 0; i < reference_stack.size(); i++)
            {
                while (!reference_stack[i].empty())
                {
                    alignments_reference[i].insert(0, 1, reference_stack[i].top());
                    reference_stack[i].pop();
                }
            }

            // pad the head for query
            for (int i = 0; i < query_stack.size(); i++)
            {
                while (alignments_query[i].length() < alignments_reference[0].length())
                {

                    alignments_query[i].insert(0, 1, '_');
                }
            }
            // pad the head for reference
            for (int i = 0; i < reference_stack.size(); i++)
            {
                while (alignments_reference[i].length() < alignments_query[0].length())
                {

                    alignments_reference[i].insert(0, 1, '_');
                }
            }

            // Merge alignments query and alignments reference into one data structure
            std::vector<string> alignments;
            for (int i = 0; i < querys.size(); i++)
            {
                alignments.push_back(alignments_query[i]);
            }
            for (int i = 0; i < references.size(); i++)
            {
                alignments.push_back(alignments_reference[i]);
            }

            return alignments;
        }

        std::vector<std::array<int, 5>> MultipleSequencesToProfileAlign(std::vector<string> seq, int len);

        template <typename T, size_t MAX_QRY_LENGTH, size_t MAX_REF_LENGTH>
        std::map<string, string> ReconstructTraceback(string query, string reference,
                                                      int query_start_idx, int reference_start_idx,
                                                      T (&tb_streams)[MAX_REF_LENGTH + MAX_QRY_LENGTH])
        {
            // If we also want to match the portion for the alignment exceeds the start index, use this
            // string alignment_query = query.substr(query_start_idx + 1, query.length());
            // string alignment_reference = reference.substr(reference_start_idx + 1, reference.length());
            string alignment_query = "";
            string alignment_reference = "";

            stack<char> query_stack;
            stack<char> reference_stack;

            // insert characters of query and reference into the stack, in sequence
            for (int i = 0; i < query_start_idx + 1; i++)
            {
                query_stack.push(query[i]);
            }
            for (int i = 0; i < reference_start_idx + 1; i++)
            {
                reference_stack.push(reference[i]);
            }

            // print traceback stream
            // printf("Traceback stream: ");
            // for (int i = 0; i < MAX_REFERENCE_LENGTH + MAX_QUERY_LENGTH; i++) {
            //     printf("%d ", tb_stream[i].to_int());
            // }

            T *curr_ptr = &tb_streams[0];
            // Insert the characters from the stack to the alignment strings to their beginning.
            // Iterating in order the tb_steram
            while (*curr_ptr != (T)AL_END_H)
            {
                // printf("curr_ptr: %d\n", curr_ptr->to_int());
                if (*curr_ptr == (T)AL_MMI_H)
                {
                    alignment_query = alignment_query.insert(0, 1, query_stack.top());
                    alignment_reference = alignment_reference.insert(0, 1, reference_stack.top());
                    query_stack.pop();
                    reference_stack.pop();
                }
                else if (*curr_ptr == (T)AL_INS_H)
                {
                    alignment_query = alignment_query.insert(0, 1, '_');
                    alignment_reference = alignment_reference.insert(0, 1, reference_stack.top());
                    reference_stack.pop();
                }
                else if (*curr_ptr == (T)AL_DEL_H)
                {
                    alignment_query = alignment_query.insert(0, 1, query_stack.top());
                    alignment_reference = alignment_reference.insert(0, 1, '_');
                    query_stack.pop();
                }
                else if (*curr_ptr == (T)AL_NULL_H)
                {
                    // Do nothing, AL_NULL Doesn't change the position
                }
                else
                {
                    printf("Alignment Output Iteartion End\n");
                }
                curr_ptr++;
            }

            // Finishing up concatenating the rest of the characters in the stack
            while (!query_stack.empty())
            {
                alignment_query = alignment_query.insert(0, 1, query_stack.top());
                query_stack.pop();
            }
            while (!reference_stack.empty())
            {
                alignment_reference = alignment_reference.insert(0, 1, reference_stack.top());
                reference_stack.pop();
            }

            // pad the head for query
            while (alignment_query.length() < alignment_reference.length())
            {
                alignment_query = alignment_query.insert(0, 1, '_');
            }
            // pad the head for reference
            while (alignment_reference.length() < alignment_query.length())
            {
                alignment_reference = alignment_reference.insert(0, 1, '_');
            }

            // return a dictionary (map) of aligned query and aligned reference
            map<string, string> alignments;
            alignments["query"] = alignment_query;
            alignments["reference"] = alignment_reference;
            return alignments;
        }

        template <typename T, size_t MAX_QRY_LENGTH, size_t MAX_REF_LENGTH>
        std::map<string, string> ReconstructTracebackSemiglobal(string query, string reference,
                                                                int query_start_idx, int reference_start_idx,
                                                                T (&tb_streams)[MAX_REF_LENGTH + MAX_QRY_LENGTH])
        {
            string alignment_reference = "";
            string alignment_query = "";
            string read_tb = "";
            T *curr_ptr = &tb_streams[0];
            if (query_start_idx == query.length() - 1)
            { // short-long
                char *qry_ptr = &query[query.length() - 1];
                char *ref_ptr = &reference[reference.length() - 1];

                for (int i = 0; i < reference.size() - reference_start_idx; i++)
                {
                    alignment_query = alignment_query.insert(0, 1, '_');
                    alignment_reference = alignment_reference.insert(0, 1, *ref_ptr--);
                }

                // cout << alignment_query << endl;
                // cout << alignment_reference << endl;
                // cout << endl;

                while (*curr_ptr != (T)AL_END_H)
                {
                    // check if ref_ptr or qry_ptr points to the first character of the string, break the loop if they do
                    if (*qry_ptr == '\0' || '\0' == *ref_ptr)
                    {
                        break;
                    }

                    // printf("curr_ptr: %d\n", curr_ptr->to_int());
                    if (*curr_ptr == (T)AL_MMI_H)
                    {
                        alignment_query = alignment_query.insert(0, 1, *qry_ptr--);
                        alignment_reference = alignment_reference.insert(0, 1, *ref_ptr--);
                    }
                    else if (*curr_ptr == (T)AL_INS_H)
                    {
                        alignment_query = alignment_query.insert(0, 1, '_');
                        alignment_reference = alignment_reference.insert(0, 1, *ref_ptr--);
                    }
                    else if (*curr_ptr == (T)AL_DEL_H)
                    {
                        alignment_query = alignment_query.insert(0, 1, *qry_ptr--);
                        alignment_reference = alignment_reference.insert(0, 1, '_');
                    }
                    else if (*curr_ptr == (T)AL_NULL_H)
                    {
                        // Do nothing, AL_NULL Doesn't change the position
                    }
                    else
                    {
                        printf("Alignment Output Iteartion End\n");
                    }
                    read_tb.append(nav_to_string(*curr_ptr));
                    curr_ptr++;

                    // cout << alignment_query << endl;
                    // cout << alignment_reference << endl;
                    // cout << read_tb << endl;
                    // cout << endl;
                }

                // finish the alignment
                while (*qry_ptr != '\0')
                {
                    alignment_query = alignment_query.insert(0, 1, *qry_ptr--);
                }
                while (*qry_ptr != '\0'){
                    alignment_reference = alignment_reference.insert(0, 1, *ref_ptr--);
                }

                while (*ref_ptr != '\0'){
                    alignment_query = alignment_query.insert(0, 1, '_');
                    alignment_reference = alignment_reference.insert(0, 1, *ref_ptr--);
                }

            }
            else if (reference_start_idx == reference.length() - 1)
            { // long-short
            }
            else
            {
                std::runtime_error("Alignment is not actually semiglobal, please check traceback start index");
            }
            cout << "Reconstructed" << endl;

            map<string, string> alignments;
            alignments["query"] = alignment_query;
            alignments["reference"] = alignment_reference;
            return alignments;
        }

        template <typename T, size_t MAX_QRY_LENGTH, size_t MAX_REF_LENGTH>
        std::map<string, string> ReconstructTracebackOverlapSuffixPrefix(string query, string reference,
                                                                int query_start_idx, int reference_start_idx,
                                                                T (&tb_streams)[MAX_REF_LENGTH + MAX_QRY_LENGTH])
        {
            string alignment_reference = "";
            string alignment_query = "";
            string read_tb = "";
            T *curr_ptr = &tb_streams[0];

            char *qry_ptr = &query[query.length() - 1];
            char *ref_ptr = &reference[reference.length() - 1];

            for (int i = 0; i < query.size() - query_start_idx; i++)
            {
                alignment_query = alignment_query.insert(0, 1, *qry_ptr--);
                alignment_reference = alignment_reference.insert(0, 1, '_');
            }





                // cout << alignment_query << endl;
                // cout << alignment_reference << endl;
                // cout << endl;

                while (*curr_ptr != (T)AL_END_H)
                {
                    // check if ref_ptr or qry_ptr points to the first character of the string, break the loop if they do
                    if (*qry_ptr == '\0' || '\0' == *ref_ptr)
                    {
                        break;
                    }

                    // printf("curr_ptr: %d\n", curr_ptr->to_int());
                    if (*curr_ptr == (T)AL_MMI_H)
                    {
                        alignment_query = alignment_query.insert(0, 1, *qry_ptr--);
                        alignment_reference = alignment_reference.insert(0, 1, *ref_ptr--);
                    }
                    else if (*curr_ptr == (T)AL_INS_H)
                    {
                        alignment_query = alignment_query.insert(0, 1, '_');
                        alignment_reference = alignment_reference.insert(0, 1, *ref_ptr--);
                    }
                    else if (*curr_ptr == (T)AL_DEL_H)
                    {
                        alignment_query = alignment_query.insert(0, 1, *qry_ptr--);
                        alignment_reference = alignment_reference.insert(0, 1, '_');
                    }
                    else if (*curr_ptr == (T)AL_NULL_H)
                    {
                        // Do nothing, AL_NULL Doesn't change the position
                    }
                    else
                    {
                        printf("Alignment Output Iteartion End\n");
                    }
                    read_tb.append(nav_to_string(*curr_ptr));
                    curr_ptr++;

                    // cout << alignment_query << endl;
                    // cout << alignment_reference << endl;
                    // cout << read_tb << endl;
                    // cout << endl;
                }

            while (*ref_ptr != '\0'){
                alignment_query = alignment_query.insert(0, 1, '_');
                alignment_reference = alignment_reference.insert(0, 1, *ref_ptr--);
            }

            cout << "Reconstructed" << endl;

            map<string, string> alignments;
            alignments["query"] = alignment_query;
            alignments["reference"] = alignment_reference;
            return alignments;
        }

        template <typename T, int LEN>
        string NavigationToString(T nav[LEN])
        {
            string str = "";
            for (int i = 0; i < LEN; i++)
            {
                str.append(HostUtils::Sequence::nav_to_string(nav[i]));
            }
            return str;
        }

        template <typename T, int NB, int MAX_QRY_LENGTH, int MAX_REF_LENGTH>
        array<map<string, string>, NB> ReconstructTracebackBlocks(string query[NB], string reference[NB],
                                                                  int query_start_idx[NB], int reference_start_idx[NB],
                                                                  T (&tb_streams)[NB][MAX_REF_LENGTH + MAX_QRY_LENGTH])
        {
            // declare the result data structure
            array<map<string, string>, NB> alignments;

            for (int i = 0; i < NB; i++)
            {
                alignments[i] = HostUtils::Sequence::ReconstructTraceback<T, MAX_QRY_LENGTH, MAX_REF_LENGTH>(query[i], reference[i],
                                                                                                             query_start_idx[i], reference_start_idx[i], tb_streams[i]);
            }

            return alignments;
        }

        template <typename T, int NB, int MAX_QRY_LENGTH, int MAX_REF_LENGTH>
        array<map<string, string>, NB> ReconstructTracebackSemiglobalBlocks(string query[NB], string reference[NB],
                                                                            int query_start_idx[NB], int reference_start_idx[NB],
                                                                            T (&tb_streams)[NB][MAX_REF_LENGTH + MAX_QRY_LENGTH])
        {

            // declare the result data structure
            array<map<string, string>, NB> alignments;

            for (int i = 0; i < NB; i++)
            {
                alignments[i] = HostUtils::Sequence::ReconstructTracebackSemiglobal<T, MAX_QRY_LENGTH, MAX_REF_LENGTH>(query[i], reference[i],
                                                                                                                       query_start_idx[i], reference_start_idx[i], tb_streams[i]);
            }

            return alignments;
        }


        template <typename T, int NB, int MAX_QRY_LENGTH, int MAX_REF_LENGTH>
        array<map<string, string>, NB> ReconstructTracebackOverlapSuffixPrefixBlocks(string query[NB], string reference[NB],
                                                                            int query_start_idx[NB], int reference_start_idx[NB],
                                                                            T (&tb_streams)[NB][MAX_REF_LENGTH + MAX_QRY_LENGTH])
        {

            // declare the result data structure
            array<map<string, string>, NB> alignments;

            for (int i = 0; i < NB; i++)
            {
                alignments[i] = HostUtils::Sequence::ReconstructTracebackOverlapSuffixPrefix<T, MAX_QRY_LENGTH, MAX_REF_LENGTH>(query[i], reference[i],
                                                                                                                       query_start_idx[i], reference_start_idx[i], tb_streams[i]);
            }

            return alignments;
        }

    }

    // template <typename T, int N_BLOCKS, typename CHAR_T, typename IDX_T, typename TBR_T>
    // std::vector<std::vector<T>>align_all_sequences(void (*seq_align_multiple_static)(const std::string&, const std::string&),
    //     std::vector<string> querys, std::vector<string> references){
    //     const size_t num_pairs;

    //     if (querys.size() != references.size()){
    //         // raise some runtime exception
    //         throw std::runtime_error("Different number of query and reference seqences. ");
    //     } else {
    //         num_pairs = querys.size();
    //     }

    //     for (int i = 0; i < num_pairs; i += N_BLOCKS){
    //         for (int j = 0; j < N_BLOCKS; j++){
    //             // align kernel once

    //         }
    //     }
    // }

    // Need ot determine the type of the character so that the traceback can be reconstructed
    // void align_all_sequences(std::vector<string> querys, std::vector<string> references, file){

    // }
}

/*
DLDDLLDDDDDUDDUUDDDDUDUDUDDUUDDDDDDDDDDDDDD******?????**?*????***?????**?????***?????*********?*????**?????**
DLDDLLDDDDDUDDUUDDDDUDUDUDDUUDDDDDDDDDDDDD

DDDLDDLDDLLLDDDDDLDLDDDUDDDLDDDDDDLDLDLDDLDDDUDDDDDD
DDDLDDLDDLLLDDDDDLDLDDDUDDDLDDDDDDLDLDLDDLDDDUDDDDDD*?**?*???****???**????****???**?******?*???**????**
*/

#endif // !HOST_H