Ticket #1: functions.h

/*
 * functions.h
 *
 *  Created on: Mar 3, 2014
 *      Author: Ismail
 */

#include <time.h>
#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <fstream>
#include <string>
#include <vector>
#include <cstring>
#include <cmath>
#include <math.h>
#include <map>
#include <stdio.h>
#include <ctype.h>
#include <sstream>
#include <stdlib.h>
#include <limits>
#include <float.h>
using namespace std;

#ifndef FUNCTIONS_H_
#define FUNCTIONS_H_

class functions {
public:
        //functions();
        //virtual ~functions();

        functions(){}

        ~functions(){}

        string x1, x2, x3;

        // performance measures
        double ppv;
        double sensitivity;
        double score_pred;
        double score_known;
        double mea;
        string pred1, pred2, pred3, known1, known2, known3;

        map<string, double> emission_count_temp;
        vector<vector<int> > seq_combinations;
        vector<int> rowvec3;
        double neg_threshold = -10000;
        double almost_zero = 10e-10;

        string prev_state;
        int dim1, dim2, dim3;
        vector<long> states;
        vector<string> seqs;
        string nucleotides = "ACGU.";
        string trio = "000";
        string trio0 = "000";
        string trio1 = "000";
        string trio2 = "000";
        string trio_seq = "000";
        map <string, double> emission_prob;
        vector<map<string, double> > emission_prob_collection;
        map <string, double> transition_from_state_count;
        map <long, int> emission_in_states_count;
        map <long, int> emission_to_states_count;

        // emission and state trans prob. read
        vector <vector<string> > state_transition;
        vector <vector<double> > state_transition_double;
        vector <vector<string> > emission;
        vector <vector<double> > emission_double;
        vector<string> rowvec4;
        vector<double> rowvec5;
        vector<string> seq_alpha;
        map<string, int> alpha_indexes_emissions;
        map<string, double> emission_indexes;
        map<string, int> indexes_states;

        string x1_aligned;
        string x2_aligned;
        string x3_aligned;
        string x1_aligned_reverse;
        string x2_aligned_reverse;
        string x3_aligned_reverse;

        char gap = '.';
        //char gap = '-';

        //vector<map<int, bool> > Table; //The vector of maps of ints and booleans
        //map<int, bool> Temp; //A temporary map to be used below
        //Table.push_back(Temp); //Add a map<int, bool> to the vector
        //Table.at(0)[1337/*key for the map*/] = true; //Value for the map

        map <string, double> rowmap;
        vector<map<string, double> > prob_trans_to_001;
        vector<map<string, double> > prob_trans_to_010;
        vector<map<string, double> > prob_trans_to_011;
        vector<map<string, double> > prob_trans_to_100;
        vector<map<string, double> > prob_trans_to_101;
        vector<map<string, double> > prob_trans_to_110;
        vector<map<string, double> > prob_trans_to_111;
        map <string, double> prev_state_to_001;         // Probability of transitions to So = 001
        map <string, double> prev_state_to_010;         // Probability of transitions to So = 010
        map <string, double> prev_state_to_011;         // Probability of transitions to So = 011
        map <string, double> prev_state_to_100;         // Probability of transitions to So = 100
        map <string, double> prev_state_to_101;         // Probability of transitions to So = 101
        map <string, double> prev_state_to_110;         // Probability of transitions to So = 110
        map <string, double> prev_state_to_111;         // Probability of transitions to So = 111

        vector<vector<double> > rowvec1;
        vector<double> rowvec2;
        vector<vector<vector<double> > > V001;
        vector<vector<vector<double> > > V010;
        vector<vector<vector<double> > > V011;
        vector<vector<vector<double> > > V100;
        vector<vector<vector<double> > > V101;
        vector<vector<vector<double> > > V110;
        vector<vector<vector<double> > > V111;

        // forward algo
        vector<vector<vector<double> > > V001_f;
        vector<vector<vector<double> > > V010_f;
        vector<vector<vector<double> > > V011_f;
        vector<vector<vector<double> > > V100_f;
        vector<vector<vector<double> > > V101_f;
        vector<vector<vector<double> > > V110_f;
        vector<vector<vector<double> > > V111_f;

        // backward algo
        vector<vector<vector<double> > > V001_b;
        vector<vector<vector<double> > > V010_b;
        vector<vector<vector<double> > > V011_b;
        vector<vector<vector<double> > > V100_b;
        vector<vector<vector<double> > > V101_b;
        vector<vector<vector<double> > > V110_b;
        vector<vector<vector<double> > > V111_b;

        //MEA
        vector<vector<vector<double> > > BM;
        string align1_reverse;
        string align2_reverse;
        string align3_reverse;
        string align1;
        string align2;
        string align3;




        //
        double max_imum(double a, double b, double c)
        {
                double maxx;
                maxx = max(a, b);
                maxx = max(maxx, c);

                return maxx;
        }

        void vector_assigner()
        {

                string number;
                stringstream strstream;
                //cout << "sizeof state_transition_double: " << state_transition_double.size() << endl;
                //cout << "sizeof state_transition_double[0]: " << state_transition_double[0].size() << endl;

                // state-transition assigner
                for (int j = 0; j < state_transition_double[0].size(); j++)
                {
                        strstream.str("");
                        strstream << states[j];
                        number = strstream.str();
                        //cout << "number: " << number << endl;
                        prev_state_to_001[number] = state_transition_double[0][j];
                        prev_state_to_010[number] = state_transition_double[1][j];
                        prev_state_to_011[number] = state_transition_double[2][j];
                        prev_state_to_100[number] = state_transition_double[3][j];
                        prev_state_to_101[number] = state_transition_double[4][j];
                        prev_state_to_110[number] = state_transition_double[5][j];
                        prev_state_to_111[number] = state_transition_double[6][j];
                        //cout << "001" << prev_state_to_001[number] << endl;
                        //cout << "010" << prev_state_to_010[number] << endl;
                        //cout << "011" << prev_state_to_011[number] << endl;
                        //cout << "100" << prev_state_to_100[number] << endl;
                        //cout << "101" << prev_state_to_101[number] << endl;
                        //cout << "110" << prev_state_to_110[number] << endl;
                        //cout << "111" << prev_state_to_111[number] << endl;

                }

                // emission probability assigner
                /*
                for (int j = 0; j < state_transition_double[0].size() - 1; j++)
                {
                        strstream.str("");
                        strstream << states[j];
                        number = strstream.str();
                        double temp1 = prev_state_to_001[number];
                        double temp2= prev_state_to_010[number];
                        double temp3= prev_state_to_011[number];
                        double temp4 = prev_state_to_100[number];
                        double temp5 = prev_state_to_101[number];
                        double temp6 = prev_state_to_110[number];
                        double temp7 = prev_state_to_111[number];
                }
                */
        }


        // reads the prob files
        int string_parser_prob(char filename[])
        {

                const int MAX_CHARS_PER_LINE = 512;
                const int MAX_TOKENS_PER_LINE = 20;
                const char* const DELIMITER = " \t,";

                int state_trans_row_index = 0;
                int emission_row_index = 0;

                // create a file-reading object
                ifstream fin;
                // fin.open("hmmparam_states.txt"); // open a file
                // fin.open("hmmparam_emissions.txt"); // open a file
                fin.open(filename); // open a file
                if (!fin.good())
                        return 1; // exit if file not found

                // read each line of the file
                while (!fin.eof())
                {
                        // read an entire line into memory
                        char buf[MAX_CHARS_PER_LINE];
                        fin.getline(buf, MAX_CHARS_PER_LINE);

                        // parse the line into blank-delimited tokens
                        int n = 0; // a for-loop index

                        // array to store memory addresses of the tokens in buf
                        const char* token[MAX_TOKENS_PER_LINE] = {}; // initialize to 0

                        // parse the line
                        token[0] = strtok(buf, DELIMITER); // first token
                        if (token[0]) // zero if line is blank
                        {
                          for (n = 1; n < MAX_TOKENS_PER_LINE; n++)
                          {
                                token[n] = strtok(0, DELIMITER); // subsequent tokens
                                if (!token[n]) break; // no more tokens
                          }
                        }

                        if (filename == "hmmparam_states_3seq.txt")
                        //if (filename == "hmmparam_states.txt")
                        {
                                //state_transition_double.push_back(rowvec2);
                                state_transition.push_back(rowvec4);
                                for (int j = 0; j < n; j++)
                                {
                                        state_transition[state_trans_row_index].push_back(token[j]);
                                        //cout << "state_transition[" << state_trans_row_index << "][" << j << "]: " << state_transition[state_trans_row_index][j];     // << endl;
                                }
                                //cout << endl;
                                state_trans_row_index++;
                        }

                        else if (filename == "hmmparam_emissions_3seq.txt")
                        //else if (filename == "hmmparam_emissions.txt")
                        {
                                emission.push_back(rowvec4);
                                for (int j = 0; j < n; j++)
                                {
                                        emission[emission_row_index].push_back(token[j]);
                                        //cout << "emission[" << emission_row_index << "][" << j << "]: " << emission[emission_row_index][j]; // << endl;
                                }
                                //cout << endl;
                                emission_row_index++;
                        }

                }

                return 0;
        }

        // finds the minimum of three values
        int minimum(int a, int b, int c)
        {
                int minm;
                minm = min(a, b);
                minm = min(minm, c);

                return minm;
        }

        double maximum(double a, double b, double c, double d, double e, double f, double g)
        {
                double maxx;
                maxx = max(a, b);
                maxx = max(maxx, c);
                maxx = max(maxx, d);
                maxx = max(maxx, e);
                maxx = max(maxx, f);
                maxx = max(maxx, g);

                return maxx;
        }

        // read text files for prob.'s
        int string_parser_emission(char filename[])
        {

                const int MAX_CHARS_PER_LINE = 512;
                const int MAX_TOKENS_PER_LINE = 5;
                const char* const DELIMITER = " \t,";

                string emission_symbol;
                int index_of_emission = 0;

                // create a file-reading object
                ifstream fin;
                // fin.open("hmmparam_states.txt"); // open a file
                // fin.open("hmmparam_emissions.txt"); // open a file
                fin.open(filename); // open a file
                if (!fin.good())
                        return 1; // exit if file not found

                // read each line of the file
                while (!fin.eof())
                {
                        // read an entire line into memory
                        char buf[MAX_CHARS_PER_LINE];
                        fin.getline(buf, MAX_CHARS_PER_LINE);

                        // parse the line into blank-delimited tokens
                        int n = 0; // a for-loop index

                        // array to store memory addresses of the tokens in buf
                        const char* token[MAX_TOKENS_PER_LINE] = {}; // initialize to 0

                        // parse the line
                        token[0] = strtok(buf, DELIMITER); // first token
                        if (token[0]) // zero if line is blank
                        {
                          for (n = 1; n < MAX_TOKENS_PER_LINE; n++)
                          {
                                token[n] = strtok(0, DELIMITER); // subsequent tokens
                                if (!token[n]) break; // no more tokens
                          }
                        }

                        //emission_symbol = token[0];
                        emission_symbol.assign(token[0], token[0] + 3);
                        emission_indexes[emission_symbol] = index_of_emission;
                        index_of_emission++;




                        //seq_combinations.push_back(rowvec3);
                        //for (int j = 0; j < n; j++)
                        //{
                        //      seq_combinations[seq_combinations_row_index].push_back(atoi(token[j]));
                        //}
                        //seq_combinations_row_index++;

                        /*
                        if (filename == "hmmparam_states.txt")
                        {
                                state_transition.push_back(rowvec);
                                for (int j = 0; j < n; j++)
                                {
                                        state_transition[state_trans_row_index].push_back(token[j]);
                                }
                                state_trans_row_index++;
                        }

                        else if (filename == "hmmparam_emissions.txt")
                        {
                                emission.push_back(rowvec);
                                for (int j = 0; j < n; j++)
                                {
                                        emission[emission_row_index].push_back(token[j]);
                                }
                                emission_row_index++;
                        }
                        */

                }

                return 0;
        }

        // read text files for prob.'s
        int string_parser(char filename[])
        {

                const int MAX_CHARS_PER_LINE = 512;
                const int MAX_TOKENS_PER_LINE = 20;
                const char* const DELIMITER = " \t,";

                int state_trans_row_index = 0;
                int emission_row_index = 0;
                int seq_combinations_row_index = 0;

                // create a file-reading object
                ifstream fin;
                // fin.open("hmmparam_states.txt"); // open a file
                // fin.open("hmmparam_emissions.txt"); // open a file
                fin.open(filename); // open a file
                if (!fin.good())
                        return 1; // exit if file not found

                // read each line of the file
                while (!fin.eof())
                {
                        // read an entire line into memory
                        char buf[MAX_CHARS_PER_LINE];
                        fin.getline(buf, MAX_CHARS_PER_LINE);

                        // parse the line into blank-delimited tokens
                        int n = 0; // a for-loop index

                        // array to store memory addresses of the tokens in buf
                        const char* token[MAX_TOKENS_PER_LINE] = {}; // initialize to 0

                        // parse the line
                        token[0] = strtok(buf, DELIMITER); // first token
                        if (token[0]) // zero if line is blank
                        {
                          for (n = 1; n < MAX_TOKENS_PER_LINE; n++)
                          {
                                token[n] = strtok(0, DELIMITER); // subsequent tokens
                                if (!token[n]) break; // no more tokens
                          }
                        }

                        seq_combinations.push_back(rowvec3);
                        for (int j = 0; j < n; j++)
                        {
                                seq_combinations[seq_combinations_row_index].push_back(atoi(token[j]));
                        }
                        seq_combinations_row_index++;

                        /*
                        if (filename == "hmmparam_states.txt")
                        {
                                state_transition.push_back(rowvec);
                                for (int j = 0; j < n; j++)
                                {
                                        state_transition[state_trans_row_index].push_back(token[j]);
                                }
                                state_trans_row_index++;
                        }

                        else if (filename == "hmmparam_emissions.txt")
                        {
                                emission.push_back(rowvec);
                                for (int j = 0; j < n; j++)
                                {
                                        emission[emission_row_index].push_back(token[j]);
                                }
                                emission_row_index++;
                        }
                        */

                }

                return 0;
        }

        // fill prob matrices
        void prob_matrix_fill()
        {
                string temp;
                //state_transition_double.push_back(rowvec2);

                //cout << "state_transition.size(): " << state_transition.size() << endl;

                // filling state transition matrix first
                for (int row = 1; row < state_transition.size(); row++)
                //for (int row = 1; row < 3; row++)
                {
                        state_transition_double.push_back(rowvec5);
                        //for (int col = 0; col < 3; col++)
                        for (int col = 0; col < state_transition[0].size(); col++)
                        {
                                temp = state_transition[row][col];
                                char *new_string = new char[temp.length() + 1];
                                strcpy(new_string, temp.c_str());
                                double var = atof(new_string);
                                state_transition_double[row-1].push_back(var);
                                //cout << "row: " << row-1 << " col: " << col << " " << state_transition_double[row-1][col] << "        ";
                        }
                        //cout << endl;
                }

                //cout << "sizeof state_transition_double: " << state_transition_double.size() << endl;


                // filling emission matrix second
                for (int row = 1; row < emission.size(); row++)
                {
                        emission_double.push_back(rowvec5);
                        for (int col = 0; col < emission[0].size(); col++)
                        {
                                temp = emission[row][col];
                                char *new_string = new char[temp.length() + 1];
                                strcpy(new_string, temp.c_str());
                                double var = atof(new_string);
                                emission_double[row-1].push_back(var);
                                //cout << emission_double[row-1][col] << " ";
                        }
                        //cout << endl;
                }
        }


        // fill in map of indexes for state transition probabilities and emission probability pairs
        void map_fill()
        {
                string temp;

                // filling map for states
                for (int i = 0; i < state_transition[0].size(); i++)
                {
                        indexes_states[state_transition[0][i]] = i + 1;
                }

                // filling map for emission
                for (int i = 1; i < emission.size(); i++)
                {
                        temp = emission[i][emission[1].size() - 1];
                        alpha_indexes_emissions[temp] = i - 1;
                }
        }

        int readfasta(char filename[])
        {
                std::ifstream input(filename);
                if(!input.good()){
                        std::cerr << "Error opening '"<<filename<<"'. Bailing out." << std::endl;
                        return -1;
                }

                int vec_index = 0;
                std::string line, name, content;
                while( std::getline( input, line ).good() )
                {
                        if( line.empty() || line[0] == '>' )     // Identifier marker
                        {
                                if( !name.empty() )
                                {
                                        // Print out what we read from the last entry
                                        // std::cout << name << " : " << content << std::endl;
                                        name.clear();
                                        seqs.push_back(content);
                                }
                                if( !line.empty() )
                                {
                                        name = line.substr(1);
                                }
                                content.clear();
                        }
                        else if( !name.empty() )
                        {
                                if( line.find(' ') != std::string::npos )               // Invalid sequence--no spaces allowed
                                {
                                        name.clear();
                                        content.clear();
                                }
                                else
                                {
                                        content += line;
                                }
                        }
                }

                if( !name.empty() )             // Print out what we read from the last entry
                {
                        // std::cout << name << " : " << content << std::endl;
                        seqs.push_back(content);
                }

                //cout << "seq1: " << seqs[0] << endl;
                //cout << "seq2: " << seqs[1] << endl;
                //cout << "seq3: " << seqs[2] << endl;

                return 0;
        }

        long dec_2_bin(long dec)
        {
                long rem,i=1,sum=0;
                do
                {
                        rem=dec%2;
                        sum=sum + (i*rem);
                        dec=dec/2;
                        i=i*10;
                }while(dec>0);
                return sum;
        }


        void standardiser()
        {

                char DUMMYCHAR = 'Y';
                for (int j = 0; j < seqs.size(); j++)
                {
                        for (int i = 0; i < seqs[j].length(); i++)
                        {
                                seqs[j][i] = toupper(seqs[j][i]);
                                if (seqs[j][i] != 'A' && seqs[j][i] != 'G' && seqs[j][i] != 'U' && seqs[j][i] != 'C' && seqs[j][i] != '.' && seqs[j][i] != '-')
                                        seqs[j][i] = 'X';
                        }
                }
        }


        void standardiser2()
        {
                char DUMMYCHAR = 'Y';

                        for (int i = 0; x1[i] != NULL; i++)
                        {
                                trio[0] = x1[i];
                                trio[1] = x2[i];
                                trio[2] = x3[i];

                                if (x1[i] == gap && x2[i] == gap && x3[i] == gap)
                                {
                                        //seqs[j].erase(seqs[j].begin()); // (str.begin()+9);

                                        x1[i] = DUMMYCHAR;
                                        x2[i] = DUMMYCHAR;
                                        x3[i] = DUMMYCHAR;
                                }
                        }

                        // you need include <algorithm> to use general algorithms like std::remove()
                        x1.erase (std::remove(x1.begin(), x1.end(), DUMMYCHAR), x1.end());
                        x2.erase (std::remove(x2.begin(), x2.end(), DUMMYCHAR), x2.end());
                        x3.erase (std::remove(x3.begin(), x3.end(), DUMMYCHAR), x3.end());
        }

        void standardiser3 ()
        {
                x1.erase (std::remove(x1.begin(), x1.end(), gap), x1.end());
                x2.erase (std::remove(x2.begin(), x2.end(), gap), x2.end());
                x3.erase (std::remove(x3.begin(), x3.end(), gap), x3.end());
        }

        void emission_prob_calc()
        {
                cout << "size of seq_comb: " << seq_combinations.size() << endl;
                string seq1, seq2, seq3;
                string number;
                stringstream strstream;
                string prev_state = "111";

                //for (int seq_count = 0; seq_count < seqs.size(); seq_count+=3)
                for (int seq_count = 0; seq_count < (seq_combinations.size() - 4); seq_count++)
                {
                        //cout << "index: " << seq_count << " " << seq_combinations[seq_count][0] << " "<< seq_combinations[seq_count][1] << " " << seq_combinations[seq_count][2] << endl;

                        //string prev_state = "111";

                        /*
                        //clearing state transition variables
                        transition_from_state_count.clear();
                        prev_state_to_001.clear();
                        prev_state_to_010.clear();
                        prev_state_to_011.clear();
                        prev_state_to_100.clear();
                        prev_state_to_101.clear();
                        prev_state_to_110.clear();
                        prev_state_to_111.clear();

                        //clearing state emission variables
                        emission_in_states_count.clear();
                        emission_count_temp.clear();
                        emission_prob.clear();
                        */

                        //seq1 = seqs[seq_count];
                        //seq2 = seqs[seq_count + 1];
                        //seq3 = seqs[seq_count + 2];

                        seq1 = seqs[seq_combinations[seq_count][0] - 1];
                        seq2 = seqs[seq_combinations[seq_count][1] - 1];
                        seq3 = seqs[seq_combinations[seq_count][2] - 1];
                        //cout << seq_combinations[seq_count][0] << endl;
                        //cout << seq_combinations[seq_count][1] << endl;
                        //cout << seq_combinations[seq_count][2] << endl;
                        //cout << "sequence1: " << seq1 << endl;
                        //cout << "sequence2: " << seq2 << endl;
                        //cout << "sequence3: " << seq3 << endl;

                        char DUMMYCHAR = 'Y';


                        // removing the triple gap state
                        for (int i = 0; i < seq1.length(); i++)
                        {
                                trio[0] = seq1[i];
                                trio[1] = seq2[i];
                                trio[2] = seq3[i];

                                if (seq1[i] == gap && seq2[i] == gap && seq3[i] == gap)
                                {
                                        //seqs[j].erase(seqs[j].begin()); // (str.begin()+9);

                                        seq1[i] = DUMMYCHAR;
                                        seq2[i] = DUMMYCHAR;
                                        seq3[i] = DUMMYCHAR;
                                        //seqs[j][i+1] = '.';
                                        //seqs[j+1][i+1] = '.';
                                        //seqs[j+2][i+1] = '.';
                                        //i--;
                                }
                        }

                        // you need include <algorithm> to use general algorithms like std::remove()
                        seq1.erase(std::remove(seq1.begin(), seq1.end(), DUMMYCHAR), seq1.end());
                        seq2.erase(std::remove(seq2.begin(), seq2.end(), DUMMYCHAR), seq2.end());
                        seq3.erase(std::remove(seq3.begin(), seq3.end(), DUMMYCHAR), seq3.end());



                        /*
                        cout << "string being processed:" << endl;
                        cout << "seq1: " << seq1 << endl;
                        cout << "seq2: " << seq2 << endl;
                        cout << "seq3: " << seq3 << endl;
                        cout << "seq1.length: " << seq1.length() << endl;
                        cout << "seq2.length: " << seq2.length() << endl;
                        cout << "seq3.length: " << seq3.length() << endl;
                        */

                        int seq_size = minimum(seq1.length(), seq2.length(), seq3.length());
                        //cout << "minimum_seq_size: " << seq_size << endl;

                        //for (int i = 0; i < seqs[seq_count].length(); i++)
                        for (int i = 0; i < seq_size; i++)
                        {



                                trio[0] = seq1[i];              //[seq_combinations[seq_count][0]][i];
                                trio[1] = seq2[i];              //[seq_combinations[seq_count][1]][i];
                                trio[2] = seq3[i];              //[seq_combinations[seq_count][2]][i];

                                if (trio == "XCC")
                                        int wait = 1;


                                //if (trio[0] != gap && trio[1] != gap && trio[2] != gap)
                                //{
                                        // calculating the counts for each emission
                                        emission_prob[trio] += 1;

                                        if (trio == "..Y")
                                        {
                                                int wait = 0;
                                        }

                                        // calculating emissions count in states
                                        if (trio[0] == gap && trio[1] == gap && trio[2] != gap)                         //INS 3
                                        {
                                                emission_in_states_count[states[0]]++;
                                                prev_state_to_001[prev_state]++;
                                                emission_count_temp[trio]++;


                                                //strstream << 111;
                                                //number = strstream.str();
                                                //cout << number << endl;
                                                //strstream.str("");
                                                //strstream << 110;
                                                //number = strstream.str();
                                                //cout << number << endl;


                                                strstream.str("");
                                                strstream << states[0];
                                                number = strstream.str();
                                                //cout << "new state: " << states[0] << endl;
                                                //cout << "number: " << states[0] << endl;
                                                prev_state = number;
                                                transition_from_state_count[number]++;
                                        }
                                        else if (trio[0] == gap && trio[1] != gap && trio[2] == gap)            //INS 2
                                        {
                                                emission_count_temp[trio]++;
                                                emission_in_states_count[states[1]]++;
                                                prev_state_to_010[prev_state]++;
                                                strstream.str("");
                                                strstream << states[1];
                                                number = strstream.str();
                                                //cout << "new state: " << states[1] << endl;
                                                //cout << "number: " << states[1] << endl;
                                                prev_state = number;
                                                transition_from_state_count[number]++;
                                        }
                                        else if (trio[0] == gap && trio[1] != gap && trio[2] != gap)            //INS 2&3
                                        {
                                                emission_count_temp[trio]++;
                                                emission_in_states_count[states[2]]++;
                                                prev_state_to_011[prev_state]++;
                                                strstream.str("");
                                                strstream << states[2];
                                                number = strstream.str();
                                                //cout << "new state: " << states[2] << endl;
                                                //cout << "number: " << states[2] << endl;
                                                prev_state = number;
                                                transition_from_state_count[number]++;
                                        }
                                        else if (trio[0] != gap && trio[1] == gap && trio[2] == gap)            //INS 1
                                        {
                                                emission_count_temp[trio]++;
                                                emission_in_states_count[states[3]]++;
                                                prev_state_to_100[prev_state]++;
                                                strstream.str("");
                                                strstream << states[3];
                                                number = strstream.str();
                                                //cout << "new state: " << states[3] << endl;
                                                //cout << "number: " << states[3] << endl;
                                                prev_state = number;
                                                transition_from_state_count[number]++;
                                        }
                                        else if (trio[0] != gap && trio[1] == gap && trio[2] != gap)            //INS 1&3
                                        {
                                                emission_count_temp[trio]++;
                                                emission_in_states_count[states[4]]++;
                                                prev_state_to_101[prev_state]++;
                                                strstream.str("");
                                                strstream << states[4];
                                                number = strstream.str();
                                                //cout << "new state: " << states[4] << endl;
                                                //cout << "number: " << states[4] << endl;
                                                prev_state = number;
                                                transition_from_state_count[number]++;
                                        }
                                        else if (trio[0] != gap && trio[1] != gap && trio[2] == gap)            //INS 1&2
                                        {
                                                emission_count_temp[trio]++;
                                                emission_in_states_count[states[5]]++;
                                                prev_state_to_110[prev_state]++;
                                                //cout << "prev_state_to_110 count: " << prev_state_to_110["110"] << endl;
                                                strstream.str("");
                                                strstream << states[5];
                                                number = strstream.str();
                                                //cout << "new state: " << states[5] << endl;
                                                //cout << "number: " << states[5] << endl;
                                                prev_state = number;
                                                transition_from_state_count[number]++;
                                        }
                                        else if (trio[0] != gap && trio[1] != gap && trio[2] != gap)            // transitions from aligned or mismatch state
                                        {
                                                emission_count_temp[trio]++;
                                                emission_in_states_count[states[6]]++;
                                                prev_state_to_111[prev_state]++;
                                                //cout << "prev_state_to_111 count: " << prev_state_to_111["110"] << endl;
                                                strstream.str("");
                                                strstream << states[6];
                                                number = strstream.str();
                                                //cout << "new state: "  << states[6] << endl;
                                                //cout << "number: " << states[6] << endl;
                                                prev_state = number;
                                                transition_from_state_count[number]++;
                                        }
                                //}
                        }


                        /*
                        for (int i = 0; i < seqs[seq_count].length(); i++)
                        {
                                trio[0] = seqs[seq_count][i];
                                trio[1] = seqs[seq_count + 1][i];
                                trio[2] = seqs[seq_count + 2][i];

                                if (trio[0] == gap && trio[1] == gap && trio[2] != gap)
                                        emission_count_temp[trio]++;
                        }
                        */


                        /*
                        cout << "110_to_001 count: " << prev_state_to_001["110"] << endl;
                        cout << "110_to_010 count: " << prev_state_to_010["110"] << endl;
                        cout << "110_to_011 count: " << prev_state_to_011["110"] << endl;
                        cout << "110_to_100 count: " << prev_state_to_100["110"] << endl;
                        cout << "110_to_101 count: " << prev_state_to_101["110"] << endl;
                        cout << "110_to_110 count: " << prev_state_to_110["110"] << endl;
                        cout << "110_to_111 count: " << prev_state_to_111["110"] << endl;
                        */

                        /*
                        // calculating probabilities
                        for( map <string, double>::const_iterator it = emission_prob.begin(); it != emission_prob.end(); ++it )
                        {
                                trio = it->first;
                                double value = it->second;

                                if (trio[0] == gap && trio[1] == gap && trio[2] != gap)                         //INS 3
                                {
                                        //cout << "emission in state 001 of with total count " << emission_in_states_count[states[0]] << endl;
                                        //cout << "emission in state 001 of: " << trio << " with count " << emission_prob[trio] << endl;
                                        emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[0]];
                                        //transitions_from_001++;
                                }
                                else if (trio[0] == gap && trio[1] != gap && trio[2] == gap)            //INS 2
                                {
                                        emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[1]];
                                        //transitions_from_010++;
                                }
                                else if (trio[0] == gap && trio[1] != gap && trio[2] != gap)            //INS 2&3
                                {
                                        emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[2]];
                                        //transitions_from_011++;
                                }
                                else if (trio[0] != gap && trio[1] == gap && trio[2] == gap)            //INS 1
                                {
                                        emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[3]];
                                        //transitions_from_100++;
                                }
                                else if (trio[0] != gap && trio[1] == gap && trio[2] != gap)            //INS 1&3
                                {
                                        emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[4]];
                                        //transitions_from_101++;
                                }
                                else if (trio[0] != gap && trio[1] != gap && trio[2] == gap)            //INS 1&2
                                {
                                        emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[5]];
                                        //transitions_from_110++;
                                }
                                else if (trio[0] != gap && trio[1] != gap && trio[2] != gap)            // transitions from aligned or mismatch state
                                {
                                        emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[6]];
                                        //transitions_from_111++;
                                }
                        }

                                        //cout << "emission prob for: " << "..U" << " in state 001: " << emission_prob["..U"] << endl;
                                        //cout << "emission prob for: " << "..A" << " in state 001: " << emission_prob["..A"] << endl;
                                        //cout << "emission prob for: " << "..C" << " in state 001: " << emission_prob["..C"] << endl;
                                        //cout << "emission prob for: " << "..G" << " in state 001: " << emission_prob["..G"] << endl;


                        //cout << "emission probability for '-C-': " << emission_prob[".C."] << endl;

                        emission_prob_collection.push_back(emission_prob);

                        //int count = 0;

                        // calculating the probabilities for state transitions
                        for(int count = 0; count < states.size(); count++)
                        {
                                strstream.str("");
                                strstream << states[count];
                                number = strstream.str();
                                //cout << "state: " << number << endl;
                                //cout << "state: " << states[count] << endl;

                                if (transition_from_state_count[number] == 0)
                                        transition_from_state_count[number] = 1;

                                //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                                //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                                prev_state_to_001[number] = prev_state_to_001[number] / transition_from_state_count[number];
                                //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                                //if (prev_state_to_001[number] < 0)
                                //      prev_state_to_001[number] = 0;

                                //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                                //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                                //cout << "prev_state_to_010[number]: " << prev_state_to_010[number] << endl;
                                prev_state_to_010[number] = prev_state_to_010[number] / transition_from_state_count[number];
                                //cout << "prev_state_to_010[number]: " << prev_state_to_010[number] << endl;
                                //if (prev_state_to_010[number] < 0)
                                //      prev_state_to_010[number] = 0;

                                //cout << "prev_state_to_011[number]: " << prev_state_to_011[number] << endl;
                                //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                                //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                                prev_state_to_011[number] = prev_state_to_011[number] / transition_from_state_count[number];
                                //cout << "prev_state_to_011[number]: " << prev_state_to_011[number] << endl;
                                //if (prev_state_to_011[number] < 0)
                                //      prev_state_to_011[number] = 0;

                                //cout << "prev_state_to_100[number]: " << prev_state_to_100[number] << endl;
                                //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                                //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                                prev_state_to_100[number] = prev_state_to_100[number] / transition_from_state_count[number];
                                //cout << "prev_state_to_100[number]: " << prev_state_to_100[number] << endl;
                                //if (prev_state_to_100[number] < 0)
                                //      prev_state_to_100[number] = 0;

                                //cout << "prev_state_to_101[number]: " << prev_state_to_101[number] << endl;
                                //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                                prev_state_to_101[number] = prev_state_to_101[number] / transition_from_state_count[number];
                                //cout << "prev_state_to_101[number]: " << prev_state_to_101[number] << endl;
                                //if (prev_state_to_101[number] < 0)
                                //      prev_state_to_101[number] = 0;

                                //cout << "prev_state_to_110[number]: " << prev_state_to_110[number] << endl;
                                //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                                prev_state_to_110[number] = prev_state_to_110[number] / transition_from_state_count[number];
                                //cout << "prev_state_to_110[number]: " << prev_state_to_110[number] << endl;
                                //if (prev_state_to_110[number] < 0)
                                //      prev_state_to_110[number] = 0;

                                //cout << "prev_state_to_111[number]: " << prev_state_to_111[number] << endl;
                                //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                                prev_state_to_111[number] = prev_state_to_111[number] / transition_from_state_count[number];
                                //cout << "prev_state_to_111[number]: " << prev_state_to_111[number] << endl;
                                //if (prev_state_to_111[number] < 0)
                                //      prev_state_to_111[number] = 0;
                        }


                        prob_trans_to_001.push_back(prev_state_to_001);
                        prob_trans_to_010.push_back(prev_state_to_010);
                        prob_trans_to_011.push_back(prev_state_to_011);
                        prob_trans_to_100.push_back(prev_state_to_100);
                        prob_trans_to_101.push_back(prev_state_to_101);
                        prob_trans_to_110.push_back(prev_state_to_110);
                        prob_trans_to_111.push_back(prev_state_to_111);


                        //cout << "prob. for transition (110->001): " << prev_state_to_001["110"] << endl;
                        //cout << "prob. for transition (110->010): " << prev_state_to_010["110"] << endl;
                        //cout << "prob. for transition (110->011): " << prev_state_to_011["110"] << endl;
                        //cout << "prob. for transition (110->100): " << prev_state_to_100["110"] << endl;
                        //cout << "prob. for transition (110->101): " << prev_state_to_101["110"] << endl;
                        //cout << "prob. for transition (110->110): " << prev_state_to_110["110"] << endl;
                        //cout << "prob. for transition (110->111): " << prev_state_to_111["110"] << endl;

                }
                */




                /*
                prob_trans_to_001.push_back(prev_state_to_001);
                prob_trans_to_010.push_back(prev_state_to_010);
                prob_trans_to_011.push_back(prev_state_to_011);
                prob_trans_to_100.push_back(prev_state_to_100);
                prob_trans_to_101.push_back(prev_state_to_101);
                prob_trans_to_110.push_back(prev_state_to_110);
                prob_trans_to_111.push_back(prev_state_to_111);
                */
                }

                cout << "end of iterations reached" << endl;


                // calculating probabilities
                for( map <string, double>::const_iterator it = emission_prob.begin(); it != emission_prob.end(); ++it )
                {
                        trio = it->first;
                        double value = it->second;

                        if (trio[0] == gap && trio[1] == gap && trio[2] != gap)                         //INS 3
                        {
                                cout << "emission in state 001 of with total count " << emission_in_states_count[states[0]] << endl;
                                cout << "emission in state 001 of: " << trio << " with count " << emission_prob[trio] << endl;
                                emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[0]];
                                //transitions_from_001++;
                        }
                        else if (trio[0] == gap && trio[1] != gap && trio[2] == gap)            //INS 2
                        {
                                emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[1]];
                                //transitions_from_010++;
                        }
                        else if (trio[0] == gap && trio[1] != gap && trio[2] != gap)            //INS 2&3
                        {
                                emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[2]];
                                //transitions_from_011++;
                        }
                        else if (trio[0] != gap && trio[1] == gap && trio[2] == gap)            //INS 1
                        {
                                emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[3]];
                                //transitions_from_100++;
                        }
                        else if (trio[0] != gap && trio[1] == gap && trio[2] != gap)            //INS 1&3
                        {
                                emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[4]];
                                //transitions_from_101++;
                        }
                        else if (trio[0] != gap && trio[1] != gap && trio[2] == gap)            //INS 1&2
                        {
                                emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[5]];
                                //transitions_from_110++;
                        }
                        else if (trio[0] != gap && trio[1] != gap && trio[2] != gap)            // transitions from aligned or mismatch state
                        {
                                emission_prob[trio] = emission_prob[trio] / emission_in_states_count[states[6]];
                                //transitions_from_111++;
                        }
                }


                cout << "emission prob for: " << "..U" << " in state 001: " << emission_prob["..U"] << endl;
                cout << "emission prob for: " << "..A" << " in state 001: " << emission_prob["..A"] << endl;
                cout << "emission prob for: " << "..C" << " in state 001: " << emission_prob["..C"] << endl;
                cout << "emission prob for: " << "..G" << " in state 001: " << emission_prob["..G"] << endl;


                //cout << "emission probability for '-C-': " << emission_prob[".C."] << endl;

                //emission_prob_collection.push_back(emission_prob);

                //int count = 0;


                // calculating the probabilities for state transitions
                for(int count = 0; count < states.size(); count++)
                {
                        strstream.str("");
                        strstream << states[count];
                        number = strstream.str();
                        //cout << "state: " << number << endl;
                        //cout << "state: " << states[count] << endl;

                        if (transition_from_state_count[number] == 0)
                                transition_from_state_count[number] = 1;

                        cout << "from state: " << number << endl;
                        cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                        cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                        prev_state_to_001[number] = prev_state_to_001[number] / transition_from_state_count[number];
                        //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                        //if (prev_state_to_001[number] < 0)
                        //      prev_state_to_001[number] = 0;

                        //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                        //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                        cout << "prev_state_to_010[number]: " << prev_state_to_010[number] << endl;
                        prev_state_to_010[number] = prev_state_to_010[number] / transition_from_state_count[number];
                        //cout << "prev_state_to_010[number]: " << prev_state_to_010[number] << endl;
                        //if (prev_state_to_010[number] < 0)
                        //      prev_state_to_010[number] = 0;

                        cout << "prev_state_to_011[number]: " << prev_state_to_011[number] << endl;
                        //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                        //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                        prev_state_to_011[number] = prev_state_to_011[number] / transition_from_state_count[number];
                        //cout << "prev_state_to_011[number]: " << prev_state_to_011[number] << endl;
                        //if (prev_state_to_011[number] < 0)
                        //      prev_state_to_011[number] = 0;

                        cout << "prev_state_to_100[number]: " << prev_state_to_100[number] << endl;
                        //cout << "prev_state_to_001[number]: " << prev_state_to_001[number] << endl;
                        //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                        prev_state_to_100[number] = prev_state_to_100[number] / transition_from_state_count[number];
                        //cout << "prev_state_to_100[number]: " << prev_state_to_100[number] << endl;
                        //if (prev_state_to_100[number] < 0)
                        //      prev_state_to_100[number] = 0;

                        cout << "prev_state_to_101[number]: " << prev_state_to_101[number] << endl;
                        //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                        prev_state_to_101[number] = prev_state_to_101[number] / transition_from_state_count[number];
                        //cout << "prev_state_to_101[number]: " << prev_state_to_101[number] << endl;
                        //if (prev_state_to_101[number] < 0)
                        //      prev_state_to_101[number] = 0;

                        cout << "prev_state_to_110[number]: " << prev_state_to_110[number] << endl;
                        //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                        prev_state_to_110[number] = prev_state_to_110[number] / transition_from_state_count[number];
                        //cout << "prev_state_to_110[number]: " << prev_state_to_110[number] << endl;
                        //if (prev_state_to_110[number] < 0)
                        //      prev_state_to_110[number] = 0;

                        cout << "prev_state_to_111[number]: " << prev_state_to_111[number] << endl;
                        //cout << "transition_from_state_count[number]: " << transition_from_state_count[number] << endl;
                        prev_state_to_111[number] = prev_state_to_111[number] / transition_from_state_count[number];
                        //cout << "prev_state_to_111[number]: " << prev_state_to_111[number] << endl;
                        //if (prev_state_to_111[number] < 0)
                        //      prev_state_to_111[number] = 0;
                }

        }

        void viterbi_matrix_initialisation()
        {
                for (int count1 = 0; count1 <= x1.length(); count1++)
                {
                        V001.push_back(rowvec1);
                        V010.push_back(rowvec1);
                        V011.push_back(rowvec1);
                        V100.push_back(rowvec1);
                        V101.push_back(rowvec1);
                        V110.push_back(rowvec1);
                        V111.push_back(rowvec1);

                        for (int count2 = 0; count2 <= x2.length(); count2++)
                        {
                                V001[count1].push_back(rowvec2);
                                V010[count1].push_back(rowvec2);
                                V011[count1].push_back(rowvec2);
                                V100[count1].push_back(rowvec2);
                                V101[count1].push_back(rowvec2);
                                V110[count1].push_back(rowvec2);
                                V111[count1].push_back(rowvec2);

                                for (int count3 = 0; count3 <= x3.length(); count3++)
                                {
                                        if (count1 == 0 || count2 == 0 || count3 == 0)
                                        {
                                                V001[count1][count2].push_back(-DBL_MAX);
                                                V010[count1][count2].push_back(-DBL_MAX);
                                                V011[count1][count2].push_back(-DBL_MAX);
                                                V100[count1][count2].push_back(-DBL_MAX);
                                                V101[count1][count2].push_back(-DBL_MAX);
                                                V110[count1][count2].push_back(-DBL_MAX);
                                                V111[count1][count2].push_back(-DBL_MAX);
                                        }
                                        else
                                        {
                                                V001[count1][count2].push_back(0);
                                                V010[count1][count2].push_back(0);
                                                V011[count1][count2].push_back(0);
                                                V100[count1][count2].push_back(0);
                                                V101[count1][count2].push_back(0);
                                                V110[count1][count2].push_back(0);
                                                V111[count1][count2].push_back(0);
                                        }
                                }
                        }
                }

                V001[0][0][0] = -DBL_MAX;
                V010[0][0][0] = -DBL_MAX;
                V011[0][0][0] = -DBL_MAX;
                V100[0][0][0] = -DBL_MAX;
                V101[0][0][0] = -DBL_MAX;
                V110[0][0][0] = -DBL_MAX;
                V111[0][0][0] = 0;

                /*cout << "printing matrices before filling: " << endl;
                cout << "printing V001 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V001[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V010 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V010[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V011 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V011[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V100 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V100[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V101 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V101[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V110 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V110[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V111 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V111[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }*/
        }

        void viterbi_3seq()
        {
                double v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7;
                double v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7;
                double v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7;
                double v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7;
                double v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7;
                double v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7;
                double v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7;

                double temp_state_111, temp_state_110, temp_state_101, temp_state_100, temp_state_011, temp_state_010, temp_state_001, tempp_state;
                double matrix_111, matrix_110, matrix_101, matrix_100, matrix_011, matrix_010, matrix_001;
                double temp_emission;
                double temp_total;


                // filling matrices
                for (dim1 = 1; dim1 <= x1.length(); dim1++)
                {
                        for (dim2 = 1; dim2 <= x2.length(); dim2++)
                        {
                                for (dim3 = 1; dim3 <= x3.length(); dim3++)
                                {

                                        // ALIGNMENT matrix
                                        trio[0] = x1[dim1-1];
                                        trio[1] = x2[dim2-1];
                                        trio[2] = x3[dim3-1];

                                        //cout << "trio: " << trio << endl;

                                        temp_state_111 = log(prev_state_to_111["111"]);
                                        temp_state_110 = log(prev_state_to_111["110"]);
                                        temp_state_101 = log(prev_state_to_111["101"]);
                                        temp_state_100 = log(prev_state_to_111["100"]);
                                        temp_state_011 = log(prev_state_to_111["11"]);
                                        temp_state_010 = log(prev_state_to_111["10"]);
                                        temp_state_001 = log(prev_state_to_111["1"]);
                                        //tempp_state = log(prev_state_to_111["001"]);

                                        matrix_111 = V111[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V110[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V101[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V100[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V011[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V010[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V001[dim1 - 1][dim2 - 1][dim3 - 1];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][6]);

                                        v111_1 = log(prev_state_to_111["111"]) + V111[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_2 = log(prev_state_to_111["110"]) + V110[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_3 = log(prev_state_to_111["101"]) + V101[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_4 = log(prev_state_to_111["100"]) + V100[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_5 = log(prev_state_to_111["11"]) + V011[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_6 = log(prev_state_to_111["10"]) + V010[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_7 = log(prev_state_to_111["1"]) + V001[dim1 - 1][dim2 - 1][dim3 - 1];

                                        // temp_total = log(emission_prob[trio]) + maximum(v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7);

                                        /*cout << "dim1: " << dim1 << " dim2: " << dim2 << " dim3: " << dim3 << endl;
                                        cout << "trio: " << trio << endl;
                                        //cout << "emission prob:" << log(emission_prob[trio]) << endl;
                                        cout << "state transition: " << log(prev_state_to_111["111"]) << endl;
                                        cout << "emission prob log: " << log(emission_double[alpha_indexes_emissions[trio]][6]) << endl;
                                        cout << "emission prob: " << emission_double[alpha_indexes_emissions[trio]][6] << endl;
                                        cout << "elements: " << v111_1 << " " << v111_2 << " "<< v111_3 << " "<< v111_4 << " "<< v111_5 << " "<< v111_6 << " "<< v111_7 << endl;*/



                                        //V111[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7);
                                        //V111[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][6]) + maximum(v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7);
                                        V111[dim1][dim2][dim3] = log(emission_double[emission_indexes[trio]][6]) + maximum(v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7);

                                        // INS3 matrix

                                        trio[0] = gap;
                                        trio[1] = gap;
                                        trio[2] = x3[dim3-1];

                                        temp_state_111 = log(prev_state_to_001["111"]);
                                        temp_state_110 = log(prev_state_to_001["110"]);
                                        temp_state_101 = log(prev_state_to_001["101"]);
                                        temp_state_100 = log(prev_state_to_001["100"]);
                                        temp_state_011 = log(prev_state_to_001["11"]);
                                        temp_state_010 = log(prev_state_to_001["10"]);
                                        temp_state_001 = log(prev_state_to_001["1"]);

                                        matrix_111 = V111[dim1][dim2][dim3 - 1];
                                        matrix_110 = V110[dim1][dim2][dim3 - 1];
                                        matrix_101 = V101[dim1][dim2][dim3 - 1];
                                        matrix_100 = V100[dim1][dim2][dim3 - 1];
                                        matrix_011 = V011[dim1][dim2][dim3 - 1];
                                        matrix_010 = V010[dim1][dim2][dim3 - 1];
                                        matrix_001 = V001[dim1][dim2][dim3 - 1];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][0]);

                                        v001_1 = log(prev_state_to_001["111"]) + V111[dim1][dim2][dim3 - 1];
                                        v001_2 = log(prev_state_to_001["110"]) + V110[dim1][dim2][dim3 - 1];
                                        v001_3 = log(prev_state_to_001["101"]) + V101[dim1][dim2][dim3 - 1];
                                        v001_4 = log(prev_state_to_001["100"]) + V100[dim1][dim2][dim3 - 1];
                                        v001_5 = log(prev_state_to_001["11"]) + V011[dim1][dim2][dim3 - 1];
                                        v001_6 = log(prev_state_to_001["10"]) + V010[dim1][dim2][dim3 - 1];
                                        v001_7 = log(prev_state_to_001["1"]) + V001[dim1][dim2][dim3 - 1];

                                        temp_total = log(emission_prob[trio]) + maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);

                                        //V001[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);
                                        //V001[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][0]) + maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);
                                        V001[dim1][dim2][dim3] = log(emission_double[emission_indexes[trio]][0]) + maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);

                                        // INS2 matrix

                                        trio[0] = gap;
                                        trio[1] = x2[dim2-1];
                                        trio[2] = gap;

                                        temp_state_111 = log(prev_state_to_010["111"]);
                                        temp_state_110 = log(prev_state_to_010["110"]);
                                        temp_state_101 = log(prev_state_to_010["101"]);
                                        temp_state_100 = log(prev_state_to_010["100"]);
                                        temp_state_011 = log(prev_state_to_010["11"]);
                                        temp_state_010 = log(prev_state_to_010["10"]);
                                        temp_state_001 = log(prev_state_to_010["1"]);

                                        matrix_111 = V111[dim1][dim2 - 1][dim3];
                                        matrix_110 = V110[dim1][dim2 - 1][dim3];
                                        matrix_101 = V101[dim1][dim2 - 1][dim3];
                                        matrix_100 = V100[dim1][dim2 - 1][dim3];
                                        matrix_011 = V011[dim1][dim2 - 1][dim3];
                                        matrix_010 = V010[dim1][dim2 - 1][dim3];
                                        matrix_001 = V001[dim1][dim2 - 1][dim3];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][1]);

                                        v010_1 = log(prev_state_to_010["111"]) + V111[dim1][dim2 - 1][dim3];
                                        v010_2 = log(prev_state_to_010["110"]) + V110[dim1][dim2 - 1][dim3];
                                        v010_3 = log(prev_state_to_010["101"]) + V101[dim1][dim2 - 1][dim3];
                                        v010_4 = log(prev_state_to_010["100"]) + V100[dim1][dim2 - 1][dim3];
                                        v010_5 = log(prev_state_to_010["11"]) + V011[dim1][dim2 - 1][dim3];
                                        v010_6 = log(prev_state_to_010["10"]) + V010[dim1][dim2 - 1][dim3];
                                        v010_7 = log(prev_state_to_010["1"]) + V001[dim1][dim2 - 1][dim3];

                                        temp_total = log(emission_prob[trio]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);

                                        //V010[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);
                                        //V010[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][1]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);
                                        V010[dim1][dim2][dim3] = log(emission_double[emission_indexes[trio]][1]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);

                                        // INS23 matrix

                                        trio[0] = gap;
                                        trio[1] = x2[dim2-1];
                                        trio[2] = x3[dim3-1];

                                        temp_state_111 = log(prev_state_to_011["111"]);
                                        temp_state_110 = log(prev_state_to_011["110"]);
                                        temp_state_101 = log(prev_state_to_011["101"]);
                                        temp_state_100 = log(prev_state_to_011["100"]);
                                        temp_state_011 = log(prev_state_to_011["11"]);
                                        temp_state_010 = log(prev_state_to_011["10"]);
                                        temp_state_001 = log(prev_state_to_011["1"]);

                                        matrix_111 = V111[dim1][dim2 - 1][dim3 - 1];
                                        matrix_110 = V110[dim1][dim2 - 1][dim3 - 1];
                                        matrix_101 = V101[dim1][dim2 - 1][dim3 - 1];
                                        matrix_100 = V100[dim1][dim2 - 1][dim3 - 1];
                                        matrix_011 = V011[dim1][dim2 - 1][dim3 - 1];
                                        matrix_010 = V010[dim1][dim2 - 1][dim3 - 1];
                                        matrix_001 = V001[dim1][dim2 - 1][dim3 - 1];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][2]);

                                        v011_1 = log(prev_state_to_011["111"]) + V111[dim1][dim2 - 1][dim3 - 1];
                                        v011_2 = log(prev_state_to_011["110"]) + V110[dim1][dim2 - 1][dim3 - 1];
                                        v011_3 = log(prev_state_to_011["101"]) + V101[dim1][dim2 - 1][dim3 - 1];
                                        v011_4 = log(prev_state_to_011["100"]) + V100[dim1][dim2 - 1][dim3 - 1];
                                        v011_5 = log(prev_state_to_011["11"]) + V011[dim1][dim2 - 1][dim3 - 1];
                                        v011_6 = log(prev_state_to_011["10"]) + V010[dim1][dim2 - 1][dim3 - 1];
                                        v011_7 = log(prev_state_to_011["1"]) + V001[dim1][dim2 - 1][dim3 - 1];

                                        temp_total = log(emission_prob[trio]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);

                                        //V011[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);
                                        //V011[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][2]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);
                                        V011[dim1][dim2][dim3] = log(emission_double[emission_indexes[trio]][2]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);

                                        // INS1 matrix

                                        trio[0] = x1[dim1-1];
                                        trio[1] = gap;
                                        trio[2] = gap;

                                        temp_state_111 = log(prev_state_to_100["111"]);
                                        temp_state_110 = log(prev_state_to_100["110"]);
                                        temp_state_101 = log(prev_state_to_100["101"]);
                                        temp_state_100 = log(prev_state_to_100["100"]);
                                        temp_state_011 = log(prev_state_to_100["11"]);
                                        temp_state_010 = log(prev_state_to_100["10"]);
                                        temp_state_001 = log(prev_state_to_100["1"]);

                                        matrix_111 = V111[dim1 - 1][dim2][dim3];
                                        matrix_110 = V110[dim1 - 1][dim2][dim3];
                                        matrix_101 = V101[dim1 - 1][dim2][dim3];
                                        matrix_100 = V100[dim1 - 1][dim2][dim3];
                                        matrix_011 = V011[dim1 - 1][dim2][dim3];
                                        matrix_010 = V010[dim1 - 1][dim2][dim3];
                                        matrix_001 = V001[dim1 - 1][dim2][dim3];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][3]);

                                        v100_1 = log(prev_state_to_100["111"]) + V111[dim1 - 1][dim2][dim3];
                                        v100_2 = log(prev_state_to_100["110"]) + V110[dim1 - 1][dim2][dim3];
                                        v100_3 = log(prev_state_to_100["101"]) + V101[dim1 - 1][dim2][dim3];
                                        v100_4 = log(prev_state_to_100["100"]) + V100[dim1 - 1][dim2][dim3];
                                        v100_5 = log(prev_state_to_100["11"]) + V011[dim1 - 1][dim2][dim3];
                                        v100_6 = log(prev_state_to_100["10"]) + V010[dim1 - 1][dim2][dim3];
                                        v100_7 = log(prev_state_to_100["1"]) + V001[dim1 - 1][dim2][dim3];

                                        temp_total = log(emission_prob[trio]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);

                                        //V100[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);
                                        //V100[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][3]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);
                                        V100[dim1][dim2][dim3] = log(emission_double[emission_indexes[trio]][3]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);

                                        // INS13 matrix

                                        trio[0] = x1[dim1-1];
                                        trio[1] = gap;
                                        trio[2] = x3[dim3-1];

                                        temp_state_111 = log(prev_state_to_101["111"]);
                                        temp_state_110 = log(prev_state_to_101["110"]);
                                        temp_state_101 = log(prev_state_to_101["101"]);
                                        temp_state_100 = log(prev_state_to_101["100"]);
                                        temp_state_011 = log(prev_state_to_101["11"]);
                                        temp_state_010 = log(prev_state_to_101["10"]);
                                        temp_state_001 = log(prev_state_to_101["1"]);

                                        matrix_111 = V111[dim1 - 1][dim2][dim3 - 1];
                                        matrix_110 = V110[dim1 - 1][dim2][dim3 - 1];
                                        matrix_101 = V101[dim1 - 1][dim2][dim3 - 1];
                                        matrix_100 = V100[dim1 - 1][dim2][dim3 - 1];
                                        matrix_011 = V011[dim1 - 1][dim2][dim3 - 1];
                                        matrix_010 = V010[dim1 - 1][dim2][dim3 - 1];
                                        matrix_001 = V001[dim1 - 1][dim2][dim3 - 1];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][4]);

                                        v101_1 = log(prev_state_to_101["111"]) + V111[dim1 - 1][dim2][dim3 - 1];
                                        v101_2 = log(prev_state_to_101["110"]) + V110[dim1 - 1][dim2][dim3 - 1];
                                        v101_3 = log(prev_state_to_101["101"]) + V101[dim1 - 1][dim2][dim3 - 1];
                                        v101_4 = log(prev_state_to_101["100"]) + V100[dim1 - 1][dim2][dim3 - 1];
                                        v101_5 = log(prev_state_to_101["11"]) + V011[dim1 - 1][dim2][dim3 - 1];
                                        v101_6 = log(prev_state_to_101["10"]) + V010[dim1 - 1][dim2][dim3 - 1];
                                        v101_7 = log(prev_state_to_101["1"]) + V001[dim1 - 1][dim2][dim3 - 1];

                                        temp_total = log(emission_prob[trio]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);

                                        //V101[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);
                                        //V101[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][4]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);
                                        V101[dim1][dim2][dim3] = log(emission_double[emission_indexes[trio]][4]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);

                                        // INS12 matrix

                                        trio[0] = x1[dim1-1];
                                        trio[1] = x2[dim2-1];
                                        trio[2] = gap;

                                        temp_state_111 = log(prev_state_to_110["111"]);
                                        temp_state_110 = log(prev_state_to_110["110"]);
                                        temp_state_101 = log(prev_state_to_110["101"]);
                                        temp_state_100 = log(prev_state_to_110["100"]);
                                        temp_state_011 = log(prev_state_to_110["11"]);
                                        temp_state_010 = log(prev_state_to_110["10"]);
                                        temp_state_001 = log(prev_state_to_110["1"]);

                                        matrix_111 = V111[dim1 - 1][dim2 - 1][dim3];
                                        matrix_110 = V110[dim1 - 1][dim2 - 1][dim3];
                                        matrix_101 = V101[dim1 - 1][dim2 - 1][dim3];
                                        matrix_100 = V100[dim1 - 1][dim2 - 1][dim3];
                                        matrix_011 = V011[dim1 - 1][dim2 - 1][dim3];
                                        matrix_010 = V010[dim1 - 1][dim2 - 1][dim3];
                                        matrix_001 = V001[dim1 - 1][dim2 - 1][dim3];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][5]);

                                        v110_1 = log(prev_state_to_110["111"]) + V111[dim1 - 1][dim2 - 1][dim3];
                                        v110_2 = log(prev_state_to_110["110"]) + V110[dim1 - 1][dim2 - 1][dim3];
                                        v110_3 = log(prev_state_to_110["101"]) + V101[dim1 - 1][dim2 - 1][dim3];
                                        v110_4 = log(prev_state_to_110["100"]) + V100[dim1 - 1][dim2 - 1][dim3];
                                        v110_5 = log(prev_state_to_110["11"]) + V011[dim1 - 1][dim2 - 1][dim3];
                                        v110_6 = log(prev_state_to_110["10"]) + V010[dim1 - 1][dim2 - 1][dim3];
                                        v110_7 = log(prev_state_to_110["1"]) + V001[dim1 - 1][dim2 - 1][dim3];

                                        temp_total = log(emission_prob[trio]) + maximum(v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7);

                                        //V110[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7);
                                        //V110[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][5]) + maximum(v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7);
                                        V110[dim1][dim2][dim3] = log(emission_double[emission_indexes[trio]][5]) + maximum(v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7);

                                }
                        }
                }

                /*
                // printing matrices
                cout << "printing matrices after filling: " << endl;
                cout << "printing V001 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V001[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V010 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V010[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V011 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V011[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V100 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V100[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V101 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V101[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V110 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V110[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "printing V111 matrix: " << endl;
                for (dim3 = 0; dim3 <= x3.length(); dim3++)
                {
                        cout << "dim3: " << dim3 << endl;
                        for (dim1 = 0; dim1 <= x1.length(); dim1++)
                        {
                                cout << "dim1: " << dim1 << endl;
                                for (dim2 = 0; dim2 <= x2.length(); dim2++)
                                {
                                        cout << V111[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }
                */


                //cout << "value: " << V100[4][3][2];


                // traceback
                dim1 = x1.length();
                dim2 = x2.length();
                dim3 = x3.length();

                prev_state = " ";
                //double tempp = maximum(V001[dim1][dim2][dim3], V010[dim1][dim2][dim3], V011[dim1][dim2][dim3], V100[dim1][dim2][dim3], V101[dim1][dim2][dim3], V110[dim1][dim2][dim3], V111[dim1][dim2][dim3]);
                //double tempp001 = V001[dim1][dim2][dim3];
                //double tempp100 = V100[dim1][dim2][dim3];
                if (V001[dim1][dim2][dim3] == maximum(V001[dim1][dim2][dim3], V010[dim1][dim2][dim3], V011[dim1][dim2][dim3], V100[dim1][dim2][dim3], V101[dim1][dim2][dim3], V110[dim1][dim2][dim3], V111[dim1][dim2][dim3]))
                        prev_state = "001";
                if (V010[dim1][dim2][dim3] == maximum(V001[dim1][dim2][dim3], V010[dim1][dim2][dim3], V011[dim1][dim2][dim3], V100[dim1][dim2][dim3], V101[dim1][dim2][dim3], V110[dim1][dim2][dim3], V111[dim1][dim2][dim3]))
                        prev_state = "010";
                if (V011[dim1][dim2][dim3] == maximum(V001[dim1][dim2][dim3], V010[dim1][dim2][dim3], V011[dim1][dim2][dim3], V100[dim1][dim2][dim3], V101[dim1][dim2][dim3], V110[dim1][dim2][dim3], V111[dim1][dim2][dim3]))
                        prev_state = "011";
                if (V100[dim1][dim2][dim3] == maximum(V001[dim1][dim2][dim3], V010[dim1][dim2][dim3], V011[dim1][dim2][dim3], V100[dim1][dim2][dim3], V101[dim1][dim2][dim3], V110[dim1][dim2][dim3], V111[dim1][dim2][dim3]))
                        prev_state = "100";
                if (V101[dim1][dim2][dim3] == maximum(V001[dim1][dim2][dim3], V010[dim1][dim2][dim3], V011[dim1][dim2][dim3], V100[dim1][dim2][dim3], V101[dim1][dim2][dim3], V110[dim1][dim2][dim3], V111[dim1][dim2][dim3]))
                        prev_state = "101";
                if (V110[dim1][dim2][dim3] == maximum(V001[dim1][dim2][dim3], V010[dim1][dim2][dim3], V011[dim1][dim2][dim3], V100[dim1][dim2][dim3], V101[dim1][dim2][dim3], V110[dim1][dim2][dim3], V111[dim1][dim2][dim3]))
                        prev_state = "110";
                if (V111[dim1][dim2][dim3] == maximum(V001[dim1][dim2][dim3], V010[dim1][dim2][dim3], V011[dim1][dim2][dim3], V100[dim1][dim2][dim3], V101[dim1][dim2][dim3], V110[dim1][dim2][dim3], V111[dim1][dim2][dim3]))
                        prev_state = "111";

                while(dim1 >= 1 && dim2 >= 1 && dim3 >= 1)
                {

                        double temp_max = 0;

                        if (prev_state == "001")        // Ins3
                        {
                                v001_1 = log(prev_state_to_001["111"]) + V111[dim1][dim2][dim3 - 1];
                                v001_2 = log(prev_state_to_001["110"]) + V110[dim1][dim2][dim3 - 1];
                                v001_3 = log(prev_state_to_001["101"]) + V101[dim1][dim2][dim3 - 1];
                                v001_4 = log(prev_state_to_001["100"]) + V100[dim1][dim2][dim3 - 1];
                                v001_5 = log(prev_state_to_001["11"]) + V011[dim1][dim2][dim3 - 1];
                                v001_6 = log(prev_state_to_001["10"]) + V010[dim1][dim2][dim3 - 1];
                                v001_7 = log(prev_state_to_001["1"]) + V001[dim1][dim2][dim3 - 1];

                                temp_max = maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);

                                if (fabs(v001_1 - temp_max) < almost_zero)
                                {
                                        prev_state = "111";
                                }
                                else if (fabs(v001_2 - temp_max) < almost_zero)
                                {
                                        prev_state = "110";
                                }
                                else if (fabs(v001_3 - temp_max) < almost_zero)
                                {
                                        prev_state = "101";
                                }
                                else if (fabs(v001_4 - temp_max) < almost_zero)
                                {
                                        prev_state = "100";
                                }
                                else if (fabs(v001_5 - temp_max) < almost_zero)
                                {
                                        prev_state = "011";
                                }
                                else if (fabs(v001_6 - temp_max) < almost_zero)
                                {
                                        prev_state = "010";
                                }
                                else if (fabs(v001_7 - temp_max) < almost_zero)
                                {
                                        prev_state = "001";
                                }
                                else
                                {
                                        cout << "No transition possible." << endl;
                                }

                                x1_aligned.push_back(gap);
                                x2_aligned.push_back(gap);
                                x3_aligned.push_back(x3[dim3-1]);
                                dim3--;
                        }
                        else if (prev_state == "010")   //Ins2
                        {
                                v010_1 = log(prev_state_to_010["111"]) + V111[dim1][dim2 - 1][dim3];
                                v010_2 = log(prev_state_to_010["110"]) + V110[dim1][dim2 - 1][dim3];
                                v010_3 = log(prev_state_to_010["101"]) + V101[dim1][dim2 - 1][dim3];
                                v010_4 = log(prev_state_to_010["100"]) + V100[dim1][dim2 - 1][dim3];
                                v010_5 = log(prev_state_to_010["11"]) + V011[dim1][dim2 - 1][dim3];
                                v010_6 = log(prev_state_to_010["10"]) + V010[dim1][dim2 - 1][dim3];
                                v010_7 = log(prev_state_to_010["1"]) + V001[dim1][dim2 - 1][dim3];

                                temp_max = maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);

                                if (fabs(v010_1 - temp_max) < almost_zero)
                                {
                                        prev_state = "111";
                                }
                                else if (fabs(v010_2 - temp_max) < almost_zero)
                                {
                                        prev_state = "110";
                                }
                                else if (fabs(v010_3 - temp_max) < almost_zero)
                                {
                                        prev_state = "101";
                                }
                                else if (fabs(v010_4 - temp_max) < almost_zero)
                                {
                                        prev_state = "100";
                                }
                                else if (fabs(v010_5 - temp_max) < almost_zero)
                                {
                                        prev_state = "011";
                                }
                                else if (fabs(v010_6 - temp_max) < almost_zero)
                                {
                                        prev_state = "010";
                                }
                                else if (fabs(v010_7 - temp_max) < almost_zero)
                                {
                                        prev_state = "001";
                                }
                                else
                                {
                                        cout << "No transition possible." << endl;
                                }

                                x1_aligned.push_back(gap);
                                x2_aligned.push_back(x2[dim2-1]);
                                x3_aligned.push_back(gap);
                                dim2--;
                        }
                        else if (prev_state == "011")   //Ins23
                        {
                                v011_1 = log(prev_state_to_011["111"]) + V111[dim1][dim2 - 1][dim3 - 1];
                                v011_2 = log(prev_state_to_011["110"]) + V110[dim1][dim2 - 1][dim3 - 1];
                                v011_3 = log(prev_state_to_011["101"]) + V101[dim1][dim2 - 1][dim3 - 1];
                                v011_4 = log(prev_state_to_011["100"]) + V100[dim1][dim2 - 1][dim3 - 1];
                                v011_5 = log(prev_state_to_011["11"]) + V011[dim1][dim2 - 1][dim3 - 1];
                                v011_6 = log(prev_state_to_011["10"]) + V010[dim1][dim2 - 1][dim3 - 1];
                                v011_7 = log(prev_state_to_011["1"]) + V001[dim1][dim2 - 1][dim3 - 1];

                                temp_max = maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);

                                if (fabs(v011_1 - temp_max) < almost_zero)
                                {
                                        prev_state = "111";
                                }
                                else if (fabs(v011_2 - temp_max) < almost_zero)
                                {
                                        prev_state = "110";
                                }
                                else if (fabs(v011_3 - temp_max) < almost_zero)
                                {
                                        prev_state = "101";
                                }
                                else if (fabs(v011_4 - temp_max) < almost_zero)
                                {
                                        prev_state = "100";
                                }
                                else if (fabs(v011_5 - temp_max) < almost_zero)
                                {
                                        prev_state = "011";
                                }
                                else if (fabs(v011_6 - temp_max) < almost_zero)
                                {
                                        prev_state = "010";
                                }
                                else if (fabs(v011_7 - temp_max) - almost_zero)
                                {
                                        prev_state = "001";
                                }
                                else
                                {
                                        cout << "No transition possible." << endl;
                                }

                                x1_aligned.push_back(gap);
                                x2_aligned.push_back(x2[dim2-1]);
                                x3_aligned.push_back(x3[dim3-1]);
                                dim2--; dim3--;
                        }
                        else if (prev_state == "100")   //Ins1
                        {
                                v100_1 = log(prev_state_to_100["111"]) + V111[dim1 - 1][dim2][dim3];
                                v100_2 = log(prev_state_to_100["110"]) + V110[dim1 - 1][dim2][dim3];
                                v100_3 = log(prev_state_to_100["101"]) + V101[dim1 - 1][dim2][dim3];
                                v100_4 = log(prev_state_to_100["100"]) + V100[dim1 - 1][dim2][dim3];
                                double tempp1 = V011[dim1 - 1][dim2][dim3];
                                double tempp2 = V010[dim1 - 1][dim2][dim3];
                                double tempp3 = V001[dim1 - 1][dim2][dim3];
                                v100_5 = log(prev_state_to_100["11"]) + V011[dim1 - 1][dim2][dim3];
                                v100_6 = log(prev_state_to_100["10"]) + V010[dim1 - 1][dim2][dim3];
                                v100_7 = log(prev_state_to_100["1"]) + V001[dim1 - 1][dim2][dim3];

                                temp_max = maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);

                                if (fabs(v100_1 - temp_max) < almost_zero)
                                {
                                        prev_state = "111";
                                }
                                else if (fabs(v100_2 - temp_max) - almost_zero)
                                {
                                        prev_state = "110";
                                }
                                else if (fabs(v100_3 - temp_max) < almost_zero)
                                {
                                        prev_state = "101";
                                }
                                else if (fabs(v100_4 - temp_max) < almost_zero)
                                {
                                        prev_state = "100";
                                }
                                else if (fabs(v100_5 - temp_max) < almost_zero)
                                {
                                        prev_state = "011";
                                }
                                else if (fabs(v100_6 - temp_max) < almost_zero)
                                {
                                        prev_state = "010";
                                }
                                else if (fabs(v100_7 - temp_max) < almost_zero)
                                {
                                        prev_state = "001";
                                }
                                else
                                {
                                        cout << "No transition possible." << endl;
                                }

                                x1_aligned.push_back(x1[dim1-1]);
                                x2_aligned.push_back(gap);
                                x3_aligned.push_back(gap);
                                dim1--;
                        }
                        else if (prev_state == "101")   //Ins13
                        {
                                v101_1 = log(prev_state_to_101["111"]) + V111[dim1 - 1][dim2][dim3 - 1];
                                v101_2 = log(prev_state_to_101["110"]) + V110[dim1 - 1][dim2][dim3 - 1];
                                v101_3 = log(prev_state_to_101["101"]) + V101[dim1 - 1][dim2][dim3 - 1];
                                v101_4 = log(prev_state_to_101["100"]) + V100[dim1 - 1][dim2][dim3 - 1];
                                v101_5 = log(prev_state_to_101["11"]) + V011[dim1 - 1][dim2][dim3 - 1];
                                v101_6 = log(prev_state_to_101["10"]) + V010[dim1 - 1][dim2][dim3 - 1];
                                v101_7 = log(prev_state_to_101["1"]) + V001[dim1 - 1][dim2][dim3 - 1];

                                temp_max = maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);

                                if (fabs(v101_1 - temp_max) < almost_zero)
                                {
                                        prev_state = "111";
                                }
                                else if (fabs(v101_2 - temp_max) < almost_zero)
                                {
                                        prev_state = "110";
                                }
                                else if (fabs(v101_3 - temp_max) < almost_zero)
                                {
                                        prev_state = "101";
                                }
                                else if (fabs(v101_4 - temp_max) < almost_zero)
                                {
                                        prev_state = "100";
                                }
                                else if (fabs(v101_5 - temp_max) < almost_zero)
                                {
                                        prev_state = "011";
                                }
                                else if (fabs(v101_6 - temp_max) < almost_zero)
                                {
                                        prev_state = "010";
                                }
                                else if (fabs(v101_7 - temp_max) < almost_zero)
                                {
                                        prev_state = "001";
                                }
                                else
                                {
                                        cout << "No transition possible." << endl;
                                }

                                x1_aligned.push_back(x1[dim1-1]);
                                x2_aligned.push_back(gap);
                                x3_aligned.push_back(x3[dim3-1]);
                                dim1--; dim3--;
                        }
                        else if (prev_state == "110")   //Ins12
                        {
                                v110_1 = log(prev_state_to_110["111"]) + V111[dim1 - 1][dim2 - 1][dim3];
                                v110_2 = log(prev_state_to_110["110"]) + V110[dim1 - 1][dim2 - 1][dim3];
                                v110_3 = log(prev_state_to_110["101"]) + V101[dim1 - 1][dim2 - 1][dim3];
                                v110_4 = log(prev_state_to_110["100"]) + V100[dim1 - 1][dim2 - 1][dim3];
                                v110_5 = log(prev_state_to_110["11"]) + V011[dim1 - 1][dim2 - 1][dim3];
                                v110_6 = log(prev_state_to_110["10"]) + V010[dim1 - 1][dim2 - 1][dim3];
                                v110_7 = log(prev_state_to_110["1"]) + V001[dim1 - 1][dim2 - 1][dim3];

                                temp_max = maximum(v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7);

                                if (fabs(v110_1 - temp_max) < almost_zero)
                                {
                                        prev_state = "111";
                                }
                                else if (fabs(v110_2 - temp_max) < almost_zero)
                                {
                                        prev_state = "110";
                                }
                                else if (fabs(v110_3 - temp_max) < almost_zero)
                                {
                                        prev_state = "101";
                                }
                                else if (fabs(v110_4 - temp_max) < almost_zero)
                                {
                                        prev_state = "100";
                                }
                                else if (fabs(v110_5 - temp_max) < almost_zero)
                                {
                                        prev_state = "011";
                                }
                                else if (fabs(v110_6 - temp_max) < almost_zero)
                                {
                                        prev_state = "010";
                                }
                                else if (fabs(v110_7 - temp_max) < almost_zero)
                                {
                                        prev_state = "001";
                                }
                                else
                                {
                                        cout << "No transition possible." << endl;
                                }

                                x1_aligned.push_back(x1[dim1-1]);
                                x2_aligned.push_back(x2[dim2-1]);
                                x3_aligned.push_back(gap);
                                dim1--; dim2--;
                        }
                        else                            //      align/mismatch
                        {
                                v111_1 = log(prev_state_to_111["111"]) + V111[dim1 - 1][dim2 - 1][dim3 - 1];
                                v111_2 = log(prev_state_to_111["110"]) + V110[dim1 - 1][dim2 - 1][dim3 - 1];
                                v111_3 = log(prev_state_to_111["101"]) + V101[dim1 - 1][dim2 - 1][dim3 - 1];
                                v111_4 = log(prev_state_to_111["100"]) + V100[dim1 - 1][dim2 - 1][dim3 - 1];
                                v111_5 = log(prev_state_to_111["11"]) + V011[dim1 - 1][dim2 - 1][dim3 - 1];
                                v111_6 = log(prev_state_to_111["10"]) + V010[dim1 - 1][dim2 - 1][dim3 - 1];
                                v111_7 = log(prev_state_to_111["1"]) + V001[dim1 - 1][dim2 - 1][dim3 - 1];

                                temp_max = maximum(v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7);

                                if (fabs(v111_1 - temp_max) < almost_zero)
                                {
                                        prev_state = "111";
                                }
                                else if (fabs(v111_2 - temp_max) < almost_zero)
                                {
                                        prev_state = "110";
                                }
                                else if (fabs(v111_3 - temp_max) < almost_zero)
                                {
                                        prev_state = "101";
                                }
                                else if (fabs(v111_4 - temp_max) < almost_zero)
                                {
                                        prev_state = "100";
                                }
                                else if (fabs(v111_5 - temp_max) < almost_zero)
                                {
                                        prev_state = "011";
                                }
                                else if (fabs(v111_6 - temp_max) < almost_zero)
                                {
                                        prev_state = "010";
                                }
                                else if (fabs(v111_7 - temp_max) < almost_zero)
                                {
                                        prev_state = "001";
                                }
                                else
                                {
                                        cout << "No transition possible." << endl;
                                }

                                x1_aligned.push_back(x1[dim1-1]);
                                x2_aligned.push_back(x2[dim2-1]);
                                x3_aligned.push_back(x3[dim3-1]);
                                dim1--; dim2--; dim3--;
                        }
                }

                for (int i = (x1_aligned.length() - 1); i >= 0; i--)
                {
                        x1_aligned_reverse.push_back(x1_aligned.at(i));
                        x2_aligned_reverse.push_back(x2_aligned.at(i));
                        x3_aligned_reverse.push_back(x3_aligned.at(i));
                }
        }

/*
        void final_probability_calc()
        {
                string temp1;
                double temp2;
                // clearing all maps
                emission_prob.clear();
                prev_state_to_001.clear();
                prev_state_to_010.clear();
                prev_state_to_011.clear();
                prev_state_to_100.clear();
                prev_state_to_101.clear();
                prev_state_to_110.clear();
                prev_state_to_111.clear();

                //map<string, int> emission_prob_counter;

                // calculating emission probs from all sequences
                for (int count = 0; count < emission_prob_collection.size(); count++)
                {
                        for ( map<string, double>::const_iterator iter =  emission_prob_collection[count].begin(); iter != emission_prob_collection[count].end(); ++iter)
                        {
                                temp1 = iter->first;
                                temp2 = iter->second;
                                if (temp2 < 0)
                                        temp2 = 0;
                                emission_prob[iter->first] += temp2;            // adding probs. of the emissions for each set of sequences
                                //emission_prob_counter[iter->first]++;         // adding the count of each set of emissions
                        }
                }

                for ( map<string, double>::const_iterator iter =  emission_prob.begin(); iter != emission_prob.end(); ++iter)
                {
                        temp1 = iter->first;
                        temp2 = iter->second;
                        emission_prob[iter->first] = iter->second / emission_prob_collection.size();            // adding probs. of the emissions for each set of sequences
                        //emission_prob[iter->first] = emission_prob[iter->first]/emission_prob_counter[iter->first];           // adding probs. of the emissions for each set of sequences
                }


                // calculating state transition probs for '->001' from all sequences
                for (int count = 0; count < prob_trans_to_001.size(); count++)
                {
                        for ( map<string, double>::const_iterator iter =  prob_trans_to_001[count].begin(); iter != prob_trans_to_001[count].end(); ++iter)
                        {
                                prev_state_to_001[iter->first] += iter->second;         // adding probs. of the emissions for each set of sequences
                                //emission_prob_counter[iter->first]++;                         // adding the count of each set of emissions
                        }

                        for ( map<string, double>::const_iterator iter =  prob_trans_to_010[count].begin(); iter != prob_trans_to_010[count].end(); ++iter)
                        {
                                prev_state_to_010[iter->first] += iter->second;         // adding probs. of the emissions for each set of sequences
                                //emission_prob_counter[iter->first]++;                         // adding the count of each set of emissions
                        }

                        for ( map<string, double>::const_iterator iter =  prob_trans_to_011[count].begin(); iter != prob_trans_to_011[count].end(); ++iter)
                        {
                                prev_state_to_011[iter->first] += iter->second;         // adding probs. of the emissions for each set of sequences
                                //emission_prob_counter[iter->first]++;                         // adding the count of each set of emissions
                        }

                        for ( map<string, double>::const_iterator iter =  prob_trans_to_100[count].begin(); iter != prob_trans_to_100[count].end(); ++iter)
                        {
                                prev_state_to_100[iter->first] += iter->second;         // adding probs. of the emissions for each set of sequences
                                //emission_prob_counter[iter->first]++;                         // adding the count of each set of emissions
                        }

                        for ( map<string, double>::const_iterator iter =  prob_trans_to_101[count].begin(); iter != prob_trans_to_101[count].end(); ++iter)
                        {
                                prev_state_to_101[iter->first] += iter->second;         // adding probs. of the emissions for each set of sequences
                                //emission_prob_counter[iter->first]++;                         // adding the count of each set of emissions
                        }

                        for ( map<string, double>::const_iterator iter =  prob_trans_to_110[count].begin(); iter != prob_trans_to_110[count].end(); ++iter)
                        {
                                prev_state_to_110[iter->first] += iter->second;         // adding probs. of the emissions for each set of sequences
                                //emission_prob_counter[iter->first]++;                         // adding the count of each set of emissions
                        }

                        for ( map<string, double>::const_iterator iter =  prob_trans_to_111[count].begin(); iter != prob_trans_to_111[count].end(); ++iter)
                        {
                                prev_state_to_111[iter->first] += iter->second;         // adding probs. of the emissions for each set of sequences
                                //emission_prob_counter[iter->first]++;                         // adding the count of each set of emissions
                        }
                }


                // dividing and finding probabilities
                for ( map<string, double>::const_iterator iter =  prev_state_to_001.begin(); iter != prev_state_to_001.end(); ++iter)
                {
                        //prev_state_to_001[iter->first] = prev_state_to_001[iter->first] / prev_state_to_001.size();           // adding probs. of the emissions for each set of sequences
                        prev_state_to_001[iter->first] = prev_state_to_001[iter->first] / prob_trans_to_001.size();
                        prev_state_to_010[iter->first] = prev_state_to_010[iter->first] / prob_trans_to_010.size();             // normalising probs. of the emissions for each set of sequences
                        prev_state_to_011[iter->first] = prev_state_to_011[iter->first] / prob_trans_to_011.size();             // normalising probs. of the emissions for each set of sequences
                        prev_state_to_100[iter->first] = prev_state_to_100[iter->first] / prob_trans_to_100.size();             // normalising probs. of the emissions for each set of sequences
                        prev_state_to_101[iter->first] = prev_state_to_101[iter->first] / prob_trans_to_101.size();             // normalising probs. of the emissions for each set of sequences
                        prev_state_to_110[iter->first] = prev_state_to_110[iter->first] / prob_trans_to_110.size();             // normalising probs. of the emissions for each set of sequences
                        prev_state_to_111[iter->first] = prev_state_to_111[iter->first] / prob_trans_to_111.size();             // normalising probs. of the emissions for each set of sequences
                }


        }
        */

        void sensitivity_and_ppv_calculator(string one_known, string two_known, string three_known)
        {

                int ppv_count = 0;
                int sensitivity_count = 0;
                for (int i = 0; i < one_known.length(); i++)
                {
                        trio1[0] = one_known[i];
                        trio1[1] = two_known[i];
                        trio1[2] = three_known[i];

                        if (i > 0)
                        {
                                trio0[0] = one_known[i-1];
                                trio0[1] = two_known[i-1];
                                trio0[2] = three_known[i-1];
                        }
                        else
                                trio0 = "XXX";

                        if (i < one_known.length() - 1)
                        {
                                trio2[0] = one_known[i+1];
                                trio2[1] = two_known[i+1];
                                trio2[2] = three_known[i+1];
                        }
                        else
                                trio2 = "XXX";

                        trio_seq[0] = x1_aligned_reverse[i];
                        trio_seq[1] = x2_aligned_reverse[i];
                        trio_seq[2] = x3_aligned_reverse[i];

                        if (trio_seq == trio0 || trio_seq == trio1 || trio_seq == trio2)
                        {
                                ppv_count++;
                                sensitivity_count++;
                        }
                }

                //int a1 = x1_aligned_reverse.length();
                //int a2 = one_known.length();
                ppv = (double)ppv_count/x1_aligned_reverse.length();
                sensitivity = (double)sensitivity_count/one_known.length();

                //cout << "ppv_count: " << ppv_count << " pred_length: " << x1_aligned_reverse.length() << endl;
                //cout << "sens_count: " << sensitivity_count << " known_length: " << one_known.length() << endl;
                //cout << "ppv: " << ppv << endl;
                //cout << "sensitivity: " << sensitivity << endl;
                ppv_count = 0;
                sensitivity_count = 0;

        }

        void alignment_scorer()
        {
                score_pred = 0;
                score_known = 0;
                int match = 2;
                int ins_del = -7;
                int mismatch = -7;
                for (int i = 0; i < x1_aligned_reverse.length(); i++)
                {
                        if ((x1_aligned_reverse[i] == x2_aligned_reverse[i]) && (x1_aligned_reverse[i]  == x3_aligned_reverse[i]) && (x2_aligned_reverse[i] == x3_aligned_reverse[i]))  // match
                        {
                                score_pred = score_pred + match;
                        }
                        else if (x1_aligned_reverse[i] == gap || x2_aligned_reverse[i] == gap || x3_aligned_reverse[i] == gap)  //
                        {
                                score_pred = score_pred + ins_del;
                        }
                        else
                        {
                                score_pred = score_pred + mismatch;
                        }
                }

                for (int i = 0; i < known1.length(); i++)
                {
                        if ((known1[i] == known2[i]) && (known1[i]  == known3[i]) && (known2[i] == known3[i]))  // match
                        {
                                score_known = score_known + match;
                        }
                        else if (x1_aligned_reverse[i] == gap || x2_aligned_reverse[i] == gap || x3_aligned_reverse[i] == gap)  //
                        {
                                score_known = score_known + ins_del;
                        }
                        else
                        {
                                score_known = score_known + mismatch;
                        }
                }
        }

        void forward_algo_3seq()
        {
                // forward algo array initializer
                for (int count1 = 0; count1 <= x1.length(); count1++)
                {
                        V001_f.push_back(rowvec1);
                        V010_f.push_back(rowvec1);
                        V011_f.push_back(rowvec1);
                        V100_f.push_back(rowvec1);
                        V101_f.push_back(rowvec1);
                        V110_f.push_back(rowvec1);
                        V111_f.push_back(rowvec1);

                        for (int count2 = 0; count2 <= x2.length(); count2++)
                        {
                                V001_f[count1].push_back(rowvec2);
                                V010_f[count1].push_back(rowvec2);
                                V011_f[count1].push_back(rowvec2);
                                V100_f[count1].push_back(rowvec2);
                                V101_f[count1].push_back(rowvec2);
                                V110_f[count1].push_back(rowvec2);
                                V111_f[count1].push_back(rowvec2);

                                for (int count3 = 0; count3 <= x3.length(); count3++)
                                {
                                        if (count1 == 0 || count2 == 0 || count3 == 0)
                                        {
                                                V001_f[count1][count2].push_back(neg_threshold);
                                                V010_f[count1][count2].push_back(neg_threshold);
                                                V011_f[count1][count2].push_back(neg_threshold);
                                                V100_f[count1][count2].push_back(neg_threshold);
                                                V101_f[count1][count2].push_back(neg_threshold);
                                                V110_f[count1][count2].push_back(neg_threshold);
                                                V111_f[count1][count2].push_back(neg_threshold);
                                        }
                                        else
                                        {
                                                V001_f[count1][count2].push_back(0);
                                                V010_f[count1][count2].push_back(0);
                                                V011_f[count1][count2].push_back(0);
                                                V100_f[count1][count2].push_back(0);
                                                V101_f[count1][count2].push_back(0);
                                                V110_f[count1][count2].push_back(0);
                                                V111_f[count1][count2].push_back(0);
                                        }
                                }
                        }
                }

                V001_f[0][0][0] = neg_threshold;
                V010_f[0][0][0] = neg_threshold;
                V011_f[0][0][0] = neg_threshold;
                V100_f[0][0][0] = neg_threshold;
                V101_f[0][0][0] = neg_threshold;
                V110_f[0][0][0] = neg_threshold;
                V111_f[0][0][0] = 0;

                float a,b,c,d,e,f,g,temp3;

                double v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7;
                double v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7;
                double v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7;
                double v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7;
                double v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7;
                double v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7;
                double v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7;

                double temp_state_111, temp_state_110, temp_state_101, temp_state_100, temp_state_011, temp_state_010, temp_state_001, tempp_state;
                double matrix_111, matrix_110, matrix_101, matrix_100, matrix_011, matrix_010, matrix_001;
                double temp_emission;
                double temp_total;


                // filling matrices
                for (dim1 = 1; dim1 <= x1.length(); dim1++)
                {
                        for (dim2 = 1; dim2 <= x2.length(); dim2++)
                        {
                                for (dim3 = 1; dim3 <= x3.length(); dim3++)
                                {

                                        // ALIGNMENT matrix
                                        trio[0] = x1[dim1-1];
                                        trio[1] = x2[dim2-1];
                                        trio[2] = x3[dim3-1];

                                        temp_state_111 = log(prev_state_to_111["111"]);
                                        temp_state_110 = log(prev_state_to_111["110"]);
                                        temp_state_101 = log(prev_state_to_111["101"]);
                                        temp_state_100 = log(prev_state_to_111["100"]);
                                        temp_state_011 = log(prev_state_to_111["11"]);
                                        temp_state_010 = log(prev_state_to_111["10"]);
                                        temp_state_001 = log(prev_state_to_111["1"]);
                                        //tempp_state = log(prev_state_to_111["001"]);


                                        matrix_111 = V111_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V110_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V101_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V100_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V011_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V010_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        matrix_111 = V001_f[dim1 - 1][dim2 - 1][dim3 - 1];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][6]);

                                        v111_1 = log(prev_state_to_111["111"]) + V111_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_2 = log(prev_state_to_111["110"]) + V110_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_3 = log(prev_state_to_111["101"]) + V101_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_4 = log(prev_state_to_111["100"]) + V100_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_5 = log(prev_state_to_111["11"]) + V011_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_6 = log(prev_state_to_111["10"]) + V010_f[dim1 - 1][dim2 - 1][dim3 - 1];
                                        v111_7 = log(prev_state_to_111["1"]) + V001_f[dim1 - 1][dim2 - 1][dim3 - 1];

                                        // a = v111_1 v111_2, v111_3, v111_4, v111_5, v111_6, v111_7;
                                        // temp3 = maximum(a, b, -DBL_MAX);
                                        // b = a + b - temp3;
                                        // a = temp3;
                                        // temp3 = a + log(1 + exp(b-a));               // sum of all state transitions

                                        a = max_imum(v111_1, v111_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                        b = v111_1 + v111_2 - a;
                                        temp3 = a + log(1 + exp(b-a));
                                        a = max_imum(v111_3, temp3, -DBL_MAX);
                                        b = v111_3 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                        a = max_imum(v111_4, temp3, -DBL_MAX);
                                        b = v111_4 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                        a = max_imum(v111_5, temp3, -DBL_MAX);
                                        b = v111_5 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                        a = max_imum(v111_6, temp3, -DBL_MAX);
                                        b = v111_6 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                        a = max_imum(v111_7, temp3, -DBL_MAX);
                                        b = v111_7 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions


                                        //temp_total = log(emission_prob[trio]) + maximum(v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7);


                                        //V111[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7);
                                        //V111_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][6]) + maximum(v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7);
                                        V111_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][6]) + temp3;

                                        // INS3 matrix

                                        trio[0] = gap;
                                        trio[1] = gap;
                                        trio[2] = x3[dim3-1];

                                        temp_state_111 = log(prev_state_to_001["111"]);
                                        temp_state_110 = log(prev_state_to_001["110"]);
                                        temp_state_101 = log(prev_state_to_001["101"]);
                                        temp_state_100 = log(prev_state_to_001["100"]);
                                        temp_state_011 = log(prev_state_to_001["11"]);
                                        temp_state_010 = log(prev_state_to_001["10"]);
                                        temp_state_001 = log(prev_state_to_001["1"]);

                                        matrix_111 = V111_f[dim1][dim2][dim3 - 1];
                                        matrix_110 = V110_f[dim1][dim2][dim3 - 1];
                                        matrix_101 = V101_f[dim1][dim2][dim3 - 1];
                                        matrix_100 = V100_f[dim1][dim2][dim3 - 1];
                                        matrix_011 = V011_f[dim1][dim2][dim3 - 1];
                                        matrix_010 = V010_f[dim1][dim2][dim3 - 1];
                                        matrix_001 = V001_f[dim1][dim2][dim3 - 1];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][0]);

                                        v001_1 = log(prev_state_to_001["111"]) + V111_f[dim1][dim2][dim3 - 1];
                                        v001_2 = log(prev_state_to_001["110"]) + V110_f[dim1][dim2][dim3 - 1];
                                        v001_3 = log(prev_state_to_001["101"]) + V101_f[dim1][dim2][dim3 - 1];
                                        v001_4 = log(prev_state_to_001["100"]) + V100_f[dim1][dim2][dim3 - 1];
                                        v001_5 = log(prev_state_to_001["11"]) + V011_f[dim1][dim2][dim3 - 1];
                                        v001_6 = log(prev_state_to_001["10"]) + V010_f[dim1][dim2][dim3 - 1];
                                        v001_7 = log(prev_state_to_001["1"]) + V001_f[dim1][dim2][dim3 - 1];


                                        a = max_imum(v001_1, v001_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                        b = v001_1 + v001_2 - a;
                                        temp3 = a + log(1 + exp(b-a));
                                        a = max_imum(v001_3, temp3, -DBL_MAX);
                                        b = v001_3 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                        a = max_imum(v001_4, temp3, -DBL_MAX);
                                        b = v001_4 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                        a = max_imum(v001_5, temp3, -DBL_MAX);
                                        b = v001_5 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                        a = max_imum(v001_6, temp3, -DBL_MAX);
                                        b = v001_6 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                        a = max_imum(v001_7, temp3, -DBL_MAX);
                                        b = v001_7 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions


                                        temp_total = log(emission_prob[trio]) + maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);

                                        //V001_f[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);
                                        V001_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][0]) + temp3;

                                        // INS2 matrix

                                        trio[0] = gap;
                                        trio[1] = x2[dim2-1];
                                        trio[2] = gap;

                                        temp_state_111 = log(prev_state_to_010["111"]);
                                        temp_state_110 = log(prev_state_to_010["110"]);
                                        temp_state_101 = log(prev_state_to_010["101"]);
                                        temp_state_100 = log(prev_state_to_010["100"]);
                                        temp_state_011 = log(prev_state_to_010["11"]);
                                        temp_state_010 = log(prev_state_to_010["10"]);
                                        temp_state_001 = log(prev_state_to_010["1"]);

                                        matrix_111 = V111_f[dim1][dim2 - 1][dim3];
                                        matrix_110 = V110_f[dim1][dim2 - 1][dim3];
                                        matrix_101 = V101_f[dim1][dim2 - 1][dim3];
                                        matrix_100 = V100_f[dim1][dim2 - 1][dim3];
                                        matrix_011 = V011_f[dim1][dim2 - 1][dim3];
                                        matrix_010 = V010_f[dim1][dim2 - 1][dim3];
                                        matrix_001 = V001_f[dim1][dim2 - 1][dim3];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][1]);

                                        v010_1 = log(prev_state_to_010["111"]) + V111_f[dim1][dim2 - 1][dim3];
                                        v010_2 = log(prev_state_to_010["110"]) + V110_f[dim1][dim2 - 1][dim3];
                                        v010_3 = log(prev_state_to_010["101"]) + V101_f[dim1][dim2 - 1][dim3];
                                        v010_4 = log(prev_state_to_010["100"]) + V100_f[dim1][dim2 - 1][dim3];
                                        v010_5 = log(prev_state_to_010["11"]) + V011_f[dim1][dim2 - 1][dim3];
                                        v010_6 = log(prev_state_to_010["10"]) + V010_f[dim1][dim2 - 1][dim3];
                                        v010_7 = log(prev_state_to_010["1"]) + V001_f[dim1][dim2 - 1][dim3];

                                        a = max_imum(v010_1, v010_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                        b = v010_1 + v010_2 - a;
                                        temp3 = a + log(1 + exp(b-a));
                                        a = max_imum(v010_3, temp3, -DBL_MAX);
                                        b = v010_3 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                        a = max_imum(v010_4, temp3, -DBL_MAX);
                                        b = v010_4 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                        a = max_imum(v010_5, temp3, -DBL_MAX);
                                        b = v010_5 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                        a = max_imum(v010_6, temp3, -DBL_MAX);
                                        b = v010_6 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                        a = max_imum(v010_7, temp3, -DBL_MAX);
                                        b = v010_7 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                        temp_total = log(emission_prob[trio]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);

                                        //V010[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);
                                        //V010_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][1]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);
                                        V010_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][1]) + temp3;

                                        // INS23 matrix

                                        trio[0] = gap;
                                        trio[1] = x2[dim2-1];
                                        trio[2] = x3[dim3-1];

                                        temp_state_111 = log(prev_state_to_011["111"]);
                                        temp_state_110 = log(prev_state_to_011["110"]);
                                        temp_state_101 = log(prev_state_to_011["101"]);
                                        temp_state_100 = log(prev_state_to_011["100"]);
                                        temp_state_011 = log(prev_state_to_011["11"]);
                                        temp_state_010 = log(prev_state_to_011["10"]);
                                        temp_state_001 = log(prev_state_to_011["1"]);

                                        matrix_111 = V111_f[dim1][dim2 - 1][dim3 - 1];
                                        matrix_110 = V110_f[dim1][dim2 - 1][dim3 - 1];
                                        matrix_101 = V101_f[dim1][dim2 - 1][dim3 - 1];
                                        matrix_100 = V100_f[dim1][dim2 - 1][dim3 - 1];
                                        matrix_011 = V011_f[dim1][dim2 - 1][dim3 - 1];
                                        matrix_010 = V010_f[dim1][dim2 - 1][dim3 - 1];
                                        matrix_001 = V001_f[dim1][dim2 - 1][dim3 - 1];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][2]);

                                        v011_1 = log(prev_state_to_011["111"]) + V111_f[dim1][dim2 - 1][dim3 - 1];
                                        v011_2 = log(prev_state_to_011["110"]) + V110_f[dim1][dim2 - 1][dim3 - 1];
                                        v011_3 = log(prev_state_to_011["101"]) + V101_f[dim1][dim2 - 1][dim3 - 1];
                                        v011_4 = log(prev_state_to_011["100"]) + V100_f[dim1][dim2 - 1][dim3 - 1];
                                        v011_5 = log(prev_state_to_011["11"]) + V011_f[dim1][dim2 - 1][dim3 - 1];
                                        v011_6 = log(prev_state_to_011["10"]) + V010_f[dim1][dim2 - 1][dim3 - 1];
                                        v011_7 = log(prev_state_to_011["1"]) + V001_f[dim1][dim2 - 1][dim3 - 1];

                                        temp_total = log(emission_prob[trio]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);

                                        a = max_imum(v011_1, v011_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                        b = v011_1 + v011_2 - a;
                                        temp3 = a + log(1 + exp(b-a));
                                        a = max_imum(v011_3, temp3, -DBL_MAX);
                                        b = v011_3 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                        a = max_imum(v011_4, temp3, -DBL_MAX);
                                        b = v011_4 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                        a = max_imum(v011_5, temp3, -DBL_MAX);
                                        b = v011_5 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                        a = max_imum(v011_6, temp3, -DBL_MAX);
                                        b = v011_6 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                        a = max_imum(v011_7, temp3, -DBL_MAX);
                                        b = v011_7 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                        //V011[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);
                                        //V011_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][2]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);
                                        V011_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][2]) + temp3;

                                        // INS1 matrix

                                        trio[0] = x1[dim1-1];
                                        trio[1] = gap;
                                        trio[2] = gap;

                                        temp_state_111 = log(prev_state_to_100["111"]);
                                        temp_state_110 = log(prev_state_to_100["110"]);
                                        temp_state_101 = log(prev_state_to_100["101"]);
                                        temp_state_100 = log(prev_state_to_100["100"]);
                                        temp_state_011 = log(prev_state_to_100["11"]);
                                        temp_state_010 = log(prev_state_to_100["10"]);
                                        temp_state_001 = log(prev_state_to_100["1"]);

                                        matrix_111 = V111_f[dim1 - 1][dim2][dim3];
                                        matrix_110 = V110_f[dim1 - 1][dim2][dim3];
                                        matrix_101 = V101_f[dim1 - 1][dim2][dim3];
                                        matrix_100 = V100_f[dim1 - 1][dim2][dim3];
                                        matrix_011 = V011_f[dim1 - 1][dim2][dim3];
                                        matrix_010 = V010_f[dim1 - 1][dim2][dim3];
                                        matrix_001 = V001_f[dim1 - 1][dim2][dim3];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][3]);

                                        v100_1 = log(prev_state_to_100["111"]) + V111_f[dim1 - 1][dim2][dim3];
                                        v100_2 = log(prev_state_to_100["110"]) + V110_f[dim1 - 1][dim2][dim3];
                                        v100_3 = log(prev_state_to_100["101"]) + V101_f[dim1 - 1][dim2][dim3];
                                        v100_4 = log(prev_state_to_100["100"]) + V100_f[dim1 - 1][dim2][dim3];
                                        v100_5 = log(prev_state_to_100["11"]) + V011_f[dim1 - 1][dim2][dim3];
                                        v100_6 = log(prev_state_to_100["10"]) + V010_f[dim1 - 1][dim2][dim3];
                                        v100_7 = log(prev_state_to_100["1"]) + V001_f[dim1 - 1][dim2][dim3];

                                        a = max_imum(v100_1, v100_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                        b = v100_1 + v100_2 - a;
                                        temp3 = a + log(1 + exp(b-a));
                                        a = max_imum(v100_3, temp3, -DBL_MAX);
                                        b = v100_3 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                        a = max_imum(v100_4, temp3, -DBL_MAX);
                                        b = v100_4 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                        a = max_imum(v100_5, temp3, -DBL_MAX);
                                        b = v100_5 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                        a = max_imum(v100_6, temp3, -DBL_MAX);
                                        b = v100_6 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                        a = max_imum(v100_7, temp3, -DBL_MAX);
                                        b = v100_7 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                        temp_total = log(emission_prob[trio]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);

                                        //V100[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);
                                        //V100_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][3]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);
                                        V100_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][3]) + temp3;

                                        // INS13 matrix

                                        trio[0] = x1[dim1-1];
                                        trio[1] = gap;
                                        trio[2] = x3[dim3-1];

                                        temp_state_111 = log(prev_state_to_101["111"]);
                                        temp_state_110 = log(prev_state_to_101["110"]);
                                        temp_state_101 = log(prev_state_to_101["101"]);
                                        temp_state_100 = log(prev_state_to_101["100"]);
                                        temp_state_011 = log(prev_state_to_101["11"]);
                                        temp_state_010 = log(prev_state_to_101["10"]);
                                        temp_state_001 = log(prev_state_to_101["1"]);

                                        matrix_111 = V111_f[dim1 - 1][dim2][dim3 - 1];
                                        matrix_110 = V110_f[dim1 - 1][dim2][dim3 - 1];
                                        matrix_101 = V101_f[dim1 - 1][dim2][dim3 - 1];
                                        matrix_100 = V100_f[dim1 - 1][dim2][dim3 - 1];
                                        matrix_011 = V011_f[dim1 - 1][dim2][dim3 - 1];
                                        matrix_010 = V010_f[dim1 - 1][dim2][dim3 - 1];
                                        matrix_001 = V001_f[dim1 - 1][dim2][dim3 - 1];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][4]);

                                        v101_1 = log(prev_state_to_101["111"]) + V111_f[dim1 - 1][dim2][dim3 - 1];
                                        v101_2 = log(prev_state_to_101["110"]) + V110_f[dim1 - 1][dim2][dim3 - 1];
                                        v101_3 = log(prev_state_to_101["101"]) + V101_f[dim1 - 1][dim2][dim3 - 1];
                                        v101_4 = log(prev_state_to_101["100"]) + V100_f[dim1 - 1][dim2][dim3 - 1];
                                        v101_5 = log(prev_state_to_101["11"]) + V011_f[dim1 - 1][dim2][dim3 - 1];
                                        v101_6 = log(prev_state_to_101["10"]) + V010_f[dim1 - 1][dim2][dim3 - 1];
                                        v101_7 = log(prev_state_to_101["1"]) + V001_f[dim1 - 1][dim2][dim3 - 1];

                                        a = max_imum(v101_1, v101_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                        b = v101_1 + v101_2 - a;
                                        temp3 = a + log(1 + exp(b-a));
                                        a = max_imum(v101_3, temp3, -DBL_MAX);
                                        b = v101_3 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                        a = max_imum(v101_4, temp3, -DBL_MAX);
                                        b = v101_4 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                        a = max_imum(v101_5, temp3, -DBL_MAX);
                                        b = v101_5 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                        a = max_imum(v101_6, temp3, -DBL_MAX);
                                        b = v101_6 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                        a = max_imum(v101_7, temp3, -DBL_MAX);
                                        b = v101_7 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                        temp_total = log(emission_prob[trio]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);

                                        //V101[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);
                                        //V101_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][4]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);
                                        V101_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][4]) + temp3;

                                        // INS12 matrix

                                        trio[0] = x1[dim1-1];
                                        trio[1] = x2[dim2-1];
                                        trio[2] = gap;

                                        temp_state_111 = log(prev_state_to_110["111"]);
                                        temp_state_110 = log(prev_state_to_110["110"]);
                                        temp_state_101 = log(prev_state_to_110["101"]);
                                        temp_state_100 = log(prev_state_to_110["100"]);
                                        temp_state_011 = log(prev_state_to_110["11"]);
                                        temp_state_010 = log(prev_state_to_110["10"]);
                                        temp_state_001 = log(prev_state_to_110["1"]);

                                        matrix_111 = V111_f[dim1 - 1][dim2 - 1][dim3];
                                        matrix_110 = V110_f[dim1 - 1][dim2 - 1][dim3];
                                        matrix_101 = V101_f[dim1 - 1][dim2 - 1][dim3];
                                        matrix_100 = V100_f[dim1 - 1][dim2 - 1][dim3];
                                        matrix_011 = V011_f[dim1 - 1][dim2 - 1][dim3];
                                        matrix_010 = V010_f[dim1 - 1][dim2 - 1][dim3];
                                        matrix_001 = V001_f[dim1 - 1][dim2 - 1][dim3];

                                        temp_emission = log(emission_prob[trio]);
                                        temp_emission = log(emission_double[alpha_indexes_emissions[trio]][5]);

                                        v110_1 = log(prev_state_to_110["111"]) + V111_f[dim1 - 1][dim2 - 1][dim3];
                                        v110_2 = log(prev_state_to_110["110"]) + V110_f[dim1 - 1][dim2 - 1][dim3];
                                        v110_3 = log(prev_state_to_110["101"]) + V101_f[dim1 - 1][dim2 - 1][dim3];
                                        v110_4 = log(prev_state_to_110["100"]) + V100_f[dim1 - 1][dim2 - 1][dim3];
                                        v110_5 = log(prev_state_to_110["11"]) + V011_f[dim1 - 1][dim2 - 1][dim3];
                                        v110_6 = log(prev_state_to_110["10"]) + V010_f[dim1 - 1][dim2 - 1][dim3];
                                        v110_7 = log(prev_state_to_110["1"]) + V001_f[dim1 - 1][dim2 - 1][dim3];

                                        a = max_imum(v110_1, v110_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                        b = v110_1 + v110_2 - a;
                                        temp3 = a + log(1 + exp(b-a));
                                        a = max_imum(v110_3, temp3, -DBL_MAX);
                                        b = v110_3 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                        a = max_imum(v110_4, temp3, -DBL_MAX);
                                        b = v110_4 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                        a = max_imum(v110_5, temp3, -DBL_MAX);
                                        b = v110_5 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                        a = max_imum(v110_6, temp3, -DBL_MAX);
                                        b = v110_6 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                        a = max_imum(v110_7, temp3, -DBL_MAX);
                                        b = v110_7 + temp3 - a;
                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                        temp_total = log(emission_prob[trio]) + maximum(v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7);

                                        //V110_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][5]) + maximum(v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7);
                                        V110_f[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][5]) + temp3;

                                }
                        }
                }

                /*//printing arrays
                cout << "V111_f: " << endl;
                for (dim1 = 1; dim1 <= x1.length(); dim1++)
                {
                        cout << "Dim1: " << dim1 << endl;
                        for (dim2 = 1; dim2 <= x2.length(); dim2++)
                        {
                                cout << "Dim2: " << dim2 << endl;
                                for (dim3 = 1; dim3 <= x3.length(); dim3++)
                                {
                                        cout << V111_f[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "V100_f: " << endl;
                for (dim1 = 1; dim1 <= x1.length(); dim1++)
                {
                        cout << "Dim1: " << dim1 << endl;
                        for (dim2 = 1; dim2 <= x2.length(); dim2++)
                        {
                                cout << "Dim2: " << dim2 << endl;
                                for (dim3 = 1; dim3 <= x3.length(); dim3++)
                                {
                                        cout << V100_f[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "V110_f: " << endl;
                for (dim1 = 1; dim1 <= x1.length(); dim1++)
                {
                        cout << "Dim1: " << dim1 << endl;
                        for (dim2 = 1; dim2 <= x2.length(); dim2++)
                        {
                                cout << "Dim2: " << dim2 << endl;
                                for (dim3 = 1; dim3 <= x3.length(); dim3++)
                                {
                                        cout << V110_f[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }*/
        }


        void backward_algo_3seq()
        {
                // backward algo array initializer
                for (int count1 = 0; count1 <= x1.length(); count1++)
                {
                        V001_b.push_back(rowvec1);
                        V010_b.push_back(rowvec1);
                        V011_b.push_back(rowvec1);
                        V100_b.push_back(rowvec1);
                        V101_b.push_back(rowvec1);
                        V110_b.push_back(rowvec1);
                        V111_b.push_back(rowvec1);

                        for (int count2 = 0; count2 <= x2.length(); count2++)
                        {
                                V001_b[count1].push_back(rowvec2);
                                V010_b[count1].push_back(rowvec2);
                                V011_b[count1].push_back(rowvec2);
                                V100_b[count1].push_back(rowvec2);
                                V101_b[count1].push_back(rowvec2);
                                V110_b[count1].push_back(rowvec2);
                                V111_b[count1].push_back(rowvec2);

                                for (int count3 = 0; count3 <= x3.length(); count3++)
                                {
                                        if (count1 == x1.length() || count2 == x2.length() || count3 == x3.length())
                                        {
                                                V001_b[count1][count2].push_back(neg_threshold);
                                                V010_b[count1][count2].push_back(neg_threshold);
                                                V011_b[count1][count2].push_back(neg_threshold);
                                                V100_b[count1][count2].push_back(neg_threshold);
                                                V101_b[count1][count2].push_back(neg_threshold);
                                                V110_b[count1][count2].push_back(neg_threshold);
                                                V111_b[count1][count2].push_back(neg_threshold);
                                        }
                                        else
                                        {
                                                V001_b[count1][count2].push_back(0);
                                                V010_b[count1][count2].push_back(0);
                                                V011_b[count1][count2].push_back(0);
                                                V100_b[count1][count2].push_back(0);
                                                V101_b[count1][count2].push_back(0);
                                                V110_b[count1][count2].push_back(0);
                                                V111_b[count1][count2].push_back(0);
                                        }
                                }
                        }
                }

                V001_b[x1.length() - 1][x2.length() - 1][x3.length() - 1] = 0;
                V010_b[x1.length() - 1][x2.length() - 1][x3.length() - 1] = 0;
                V011_b[x1.length() - 1][x2.length() - 1][x3.length() - 1] = 0;
                V100_b[x1.length() - 1][x2.length() - 1][x3.length() - 1] = 0;
                V101_b[x1.length() - 1][x2.length() - 1][x3.length() - 1] = 0;
                V110_b[x1.length() - 1][x2.length() - 1][x3.length() - 1] = 0;
                V111_b[x1.length() - 1][x2.length() - 1][x3.length() - 1] = 0;

                float a,b,c,d,e,f,g,temp3;

                double v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7;
                double v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7;
                double v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7;
                double v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7;
                double v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7;
                double v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7;
                double v111_1, v111_2, v111_3, v111_4, v111_5, v111_6, v111_7;

                double emm_111, emm_110, emm_101, emm_100, emm_011, emm_010, emm_001;
                double temp_state_111, temp_state_110, temp_state_101, temp_state_100, temp_state_011, temp_state_010, temp_state_001, tempp_state;
                double matrix_111, matrix_110, matrix_101, matrix_100, matrix_011, matrix_010, matrix_001;
                double temp_emission1, temp_emission2, temp_emission3, temp_emission4, temp_emission5, temp_emission6, temp_emission7;
                double temp_total;

                char triplet_111[4] = {0};
                char triplet_110[4] = {0};
                char triplet_101[4] = {0};
                char triplet_100[4] = {0};
                char triplet_011[4] = {0};
                char triplet_010[4] = {0};
                char triplet_001[4] = {0};

                // filling matrices
                for (dim1 = x1.length()-1; dim1 >= 0; dim1--)
                {
                        for (dim2 = x2.length()-1; dim2 >= 0; dim2--)
                        {
                                for (dim3 = x3.length()-1; dim3 >= 0; dim3--)
                                {
                                                if (dim3 < (x3.length() - 1))   // 001
                                                {
                                                        triplet_001[0] = gap;
                                                        triplet_001[1] = gap;
                                                        triplet_001[2] = x3[dim3+1];
                                                        emm_001 = log(emission_double[alpha_indexes_emissions[triplet_001]][0]);
                                                }
                                                else
                                                {
                                                        triplet_001[0] = gap;
                                                        triplet_001[1] = gap;
                                                        triplet_001[2] = 'A';
                                                        emm_001 = neg_threshold;
                                                }

                                                if (dim2 < (x2.length() - 1))   // 010
                                                {
                                                        triplet_010[0] = gap;
                                                        triplet_010[1] = x2[dim2+1];
                                                        triplet_010[2] = gap;
                                                        emm_010 = log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                }
                                                else
                                                {
                                                        triplet_010[0] = gap;
                                                        triplet_010[1] = 'A';
                                                        triplet_010[2] = gap;
                                                        emm_010 = neg_threshold;
                                                }

                                                if (dim2 < (x2.length() - 1) && dim3 < (x3.length() - 1))       // 011
                                                {
                                                        triplet_011[0] = gap;
                                                        triplet_011[1] = x2[dim2+1];
                                                        triplet_011[2] = x3[dim3+1];
                                                        emm_011 = log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                }
                                                else
                                                {
                                                        triplet_011[0] = gap;
                                                        triplet_011[1] = 'A';
                                                        triplet_011[2] = 'A';
                                                        emm_011 = neg_threshold;
                                                }

                                                if (dim1 < (x1.length() - 1))   // 100
                                                {
                                                        triplet_100[0] = x1[dim1+1];
                                                        triplet_100[1] = gap;
                                                        triplet_100[2] = gap;
                                                        emm_100 = log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                }
                                                else
                                                {
                                                        triplet_100[0] = 'A';
                                                        triplet_100[1] = gap;
                                                        triplet_100[2] = gap;
                                                        emm_100 = neg_threshold;
                                                }

                                                if ((dim1 < (x1.length() - 1)) && (dim3 < (x3.length() - 1)))   // 101
                                                {
                                                        triplet_101[0] = x1[dim1+1];
                                                        triplet_101[1] = gap;
                                                        triplet_101[2] = x3[dim3+1];
                                                        emm_101 = log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                }
                                                else
                                                {
                                                        triplet_101[0] = 'A';
                                                        triplet_101[1] = gap;
                                                        triplet_101[2] = 'A';
                                                        emm_101 = neg_threshold;
                                                }

                                                if ((dim1 < (x1.length() - 1)) && (dim2 < (x2.length() - 1)))   // 110
                                                {
                                                        triplet_110[0] = x1[dim1+1];
                                                        triplet_110[1] = x2[dim2+1];
                                                        triplet_110[2] = gap;
                                                        emm_110 = log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                }
                                                else
                                                {
                                                        triplet_110[0] = 'A';
                                                        triplet_110[1] = 'A';
                                                        triplet_110[2] = gap;
                                                        emm_110 = neg_threshold;
                                                }

                                                if (dim1 < (x1.length() - 1) && dim2 < (x2.length() - 1) && dim3 < (x3.length() - 1))   // 111
                                                {
                                                        triplet_111[0] = x1[dim1+1];
                                                        triplet_111[1] = x2[dim2+1];
                                                        triplet_111[2] = x3[dim3+1];
                                                        emm_111 = log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                }
                                                else
                                                {
                                                        triplet_111[0] = 'A';
                                                        triplet_111[1] = 'A';
                                                        triplet_111[2] = 'A';
                                                        emm_111 = neg_threshold;
                                                }

                                                if (!((dim1 == x1.length()-1) && (dim2 == x2.length()-1) && (dim3 == x3.length()-1)))
                                                {
                                                        // aligned
                                                        temp_state_111 = log(prev_state_to_111["111"]);
                                                        temp_state_110 = log(prev_state_to_110["111"]);
                                                        temp_state_101 = log(prev_state_to_101["111"]);
                                                        temp_state_100 = log(prev_state_to_100["111"]);
                                                        temp_state_011 = log(prev_state_to_011["111"]);
                                                        temp_state_010 = log(prev_state_to_010["111"]);
                                                        temp_state_001 = log(prev_state_to_001["111"]);

                                                        temp_emission1 = log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                        temp_emission2 = log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                        temp_emission3 = log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                        temp_emission4 = log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                        temp_emission5 = log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                        temp_emission6 = log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                        temp_emission7 = log(emission_double[alpha_indexes_emissions[triplet_001]][0]);

                                                        matrix_111 = V111_b[dim1 + 1][dim2 + 1][dim3 + 1];
                                                        matrix_110 = V110_b[dim1 + 1][dim2 + 1][dim3];
                                                        matrix_101 = V101_b[dim1 + 1][dim2][dim3 + 1];
                                                        matrix_100 = V100_b[dim1 + 1][dim2][dim3];
                                                        matrix_011 = V011_b[dim1][dim2 + 1][dim3 + 1];
                                                        matrix_010 = V010_b[dim1][dim2 + 1][dim3];
                                                        matrix_001 = V001_b[dim1][dim2][dim3 + 1];

                                                        temp_emission1 = log(emission_prob[trio]);
                                                        //temp_emission = log(emission_double[alpha_indexes_emissions[trio]][6]);

                                                        v111_1 = log(prev_state_to_111["111"]) + V111_b[dim1 + 1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                        v111_2 = log(prev_state_to_110["111"]) + V110_b[dim1 + 1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                        v111_3 = log(prev_state_to_101["111"]) + V101_b[dim1 + 1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                        v111_4 = log(prev_state_to_100["111"]) + V100_b[dim1 + 1][dim2][dim3] + log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                        v111_5 = log(prev_state_to_011["111"]) + V011_b[dim1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                        v111_6 = log(prev_state_to_010["111"]) + V010_b[dim1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                        v111_7 = log(prev_state_to_001["111"]) + V001_b[dim1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_001]][0]);

                                                        a = max_imum(v111_1, v111_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                                        b = v111_1 + v111_2 - a;
                                                        temp3 = a + log(1 + exp(b-a));
                                                        a = max_imum(v111_3, temp3, -DBL_MAX);
                                                        b = v111_3 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                                        a = max_imum(v111_4, temp3, -DBL_MAX);
                                                        b = v111_4 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                                        a = max_imum(v111_5, temp3, -DBL_MAX);
                                                        b = v111_5 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                                        a = max_imum(v111_6, temp3, -DBL_MAX);
                                                        b = v111_6 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                                        a = max_imum(v111_7, temp3, -DBL_MAX);
                                                        b = v111_7 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                                        V111_b[dim1][dim2][dim3] = temp3;

                                                        // INS3 matrix
                                                        trio[0] = gap;
                                                        trio[1] = gap;
                                                        trio[2] = x3[dim3-1];

                                                        temp_state_111 = log(prev_state_to_111["1"]);
                                                        temp_state_110 = log(prev_state_to_110["1"]);
                                                        temp_state_101 = log(prev_state_to_101["1"]);
                                                        temp_state_100 = log(prev_state_to_100["1"]);
                                                        temp_state_011 = log(prev_state_to_011["1"]);
                                                        temp_state_010 = log(prev_state_to_010["1"]);
                                                        temp_state_001 = log(prev_state_to_001["1"]);

                                                        temp_emission1 = log(emission_prob[trio]);
                                                        temp_emission1 = log(emission_double[alpha_indexes_emissions[trio]][0]);

                                                        v001_1 = log(prev_state_to_111["1"]) + V111_b[dim1 + 1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                        v001_2 = log(prev_state_to_110["1"]) + V110_b[dim1 + 1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                        v001_3 = log(prev_state_to_101["1"]) + V101_b[dim1 + 1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                        v001_4 = log(prev_state_to_100["1"]) + V100_b[dim1 + 1][dim2][dim3] + log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                        v001_5 = log(prev_state_to_011["1"]) + V011_b[dim1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                        v001_6 = log(prev_state_to_010["1"]) + V010_b[dim1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                        v001_7 = log(prev_state_to_001["1"]) + V001_b[dim1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_001]][0]);

                                                        a = max_imum(v001_1, v001_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                                        b = v001_1 + v001_2 - a;
                                                        temp3 = a + log(1 + exp(b-a));
                                                        a = max_imum(v001_3, temp3, -DBL_MAX);
                                                        b = v001_3 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                                        a = max_imum(v001_4, temp3, -DBL_MAX);
                                                        b = v001_4 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                                        a = max_imum(v001_5, temp3, -DBL_MAX);
                                                        b = v001_5 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                                        a = max_imum(v001_6, temp3, -DBL_MAX);
                                                        b = v001_6 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                                        a = max_imum(v001_7, temp3, -DBL_MAX);
                                                        b = v001_7 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                                        temp_total = log(emission_prob[trio]) + maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);

                                                        //V001[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v001_1, v001_2, v001_3, v001_4, v001_5, v001_6, v001_7);
                                                        V001_b[dim1][dim2][dim3] = temp3;

                                                        // INS2 matrix

                                                        //trio[0] = gap;
                                                        //trio[1] = x2[dim2-1];
                                                        //trio[2] = gap;

                                                        temp_state_111 = log(prev_state_to_010["111"]);
                                                        temp_state_110 = log(prev_state_to_010["110"]);
                                                        temp_state_101 = log(prev_state_to_010["101"]);
                                                        temp_state_100 = log(prev_state_to_010["100"]);
                                                        temp_state_011 = log(prev_state_to_010["11"]);
                                                        temp_state_010 = log(prev_state_to_010["10"]);
                                                        temp_state_001 = log(prev_state_to_010["1"]);

                                                        //temp_emission1 = log(emission_prob[trio]);
                                                        //temp_emission1 = log(emission_double[alpha_indexes_emissions[trio]][1]);

                                                        v010_1 = log(prev_state_to_111["10"]) + V111_b[dim1 + 1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                        v010_2 = log(prev_state_to_110["10"]) + V110_b[dim1 + 1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                        v010_3 = log(prev_state_to_101["10"]) + V101_b[dim1 + 1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                        v010_4 = log(prev_state_to_100["10"]) + V100_b[dim1 + 1][dim2][dim3] + log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                        v010_5 = log(prev_state_to_011["10"]) + V011_b[dim1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                        v010_6 = log(prev_state_to_010["10"]) + V010_b[dim1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                        v010_7 = log(prev_state_to_001["10"]) + V001_b[dim1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_001]][0]);

                                                        a = max_imum(v010_1, v010_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                                        b = v010_1 + v010_2 - a;
                                                        temp3 = a + log(1 + exp(b-a));
                                                        a = max_imum(v010_3, temp3, -DBL_MAX);
                                                        b = v010_3 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                                        a = max_imum(v010_4, temp3, -DBL_MAX);
                                                        b = v010_4 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                                        a = max_imum(v010_5, temp3, -DBL_MAX);
                                                        b = v010_5 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                                        a = max_imum(v010_6, temp3, -DBL_MAX);
                                                        b = v010_6 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                                        a = max_imum(v010_7, temp3, -DBL_MAX);
                                                        b = v010_7 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                                        temp_total = log(emission_prob[trio]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);

                                                        //V010[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);
                                                        //V010[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][1]) + maximum(v010_1, v010_2, v010_3, v010_4, v010_5, v010_6, v010_7);
                                                        V010_b[dim1][dim2][dim3] = temp3;

                                                        // INS23 matrix

                                                        trio[0] = gap;
                                                        trio[1] = x2[dim2-1];
                                                        trio[2] = x3[dim3-1];

                                                        temp_state_111 = log(prev_state_to_011["111"]);
                                                        temp_state_110 = log(prev_state_to_011["110"]);
                                                        temp_state_101 = log(prev_state_to_011["101"]);
                                                        temp_state_100 = log(prev_state_to_011["100"]);
                                                        temp_state_011 = log(prev_state_to_011["11"]);
                                                        temp_state_010 = log(prev_state_to_011["10"]);
                                                        temp_state_001 = log(prev_state_to_011["1"]);

                                                        //temp_emission1 = log(emission_prob[trio]);
                                                        //temp_emission1 = log(emission_double[alpha_indexes_emissions[trio]][2]);

                                                        v011_1 = log(prev_state_to_111["11"]) + V111_b[dim1 + 1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                        v011_2 = log(prev_state_to_110["11"]) + V110_b[dim1 + 1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                        v011_3 = log(prev_state_to_101["11"]) + V101_b[dim1 + 1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                        v011_4 = log(prev_state_to_100["11"]) + V100_b[dim1 + 1][dim2][dim3] + log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                        v011_5 = log(prev_state_to_011["11"]) + V011_b[dim1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                        v011_6 = log(prev_state_to_010["11"]) + V010_b[dim1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                        v011_7 = log(prev_state_to_001["11"]) + V001_b[dim1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_001]][0]);

                                                        //v011_1 = log(prev_state_to_111["011"]) + V111[dim1][dim2 - 1][dim3 - 1];
                                                        //v011_2 = log(prev_state_to_111["011"]) + V110[dim1][dim2 - 1][dim3 - 1];
                                                        //v011_3 = log(prev_state_to_011["011"]) + V101[dim1][dim2 - 1][dim3 - 1];
                                                        //v011_4 = log(prev_state_to_011["011"]) + V100[dim1][dim2 - 1][dim3 - 1];
                                                        //v011_5 = log(prev_state_to_011["011"]) + V011[dim1][dim2 - 1][dim3 - 1];
                                                        //v011_6 = log(prev_state_to_011["011"]) + V010[dim1][dim2 - 1][dim3 - 1];
                                                        //v011_7 = log(prev_state_to_011["011"]) + V001[dim1][dim2 - 1][dim3 - 1];

                                                        //temp_total = log(emission_prob[trio]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);

                                                        a = max_imum(v011_1, v010_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                                        b = v011_1 + v011_2 - a;
                                                        temp3 = a + log(1 + exp(b-a));
                                                        a = max_imum(v011_3, temp3, -DBL_MAX);
                                                        b = v011_3 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                                        a = max_imum(v011_4, temp3, -DBL_MAX);
                                                        b = v011_4 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                                        a = max_imum(v011_5, temp3, -DBL_MAX);
                                                        b = v011_5 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                                        a = max_imum(v011_6, temp3, -DBL_MAX);
                                                        b = v011_6 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                                        a = max_imum(v011_7, temp3, -DBL_MAX);
                                                        b = v011_7 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                                        //V011[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);
                                                        //V011[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][2]) + maximum(v011_1, v011_2, v011_3, v011_4, v011_5, v011_6, v011_7);
                                                        //V011[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][2]) + temp3;
                                                        V011_b[dim1][dim2][dim3] = temp3;

                                                        // INS1 matrix

                                                        trio[0] = x1[dim1-1];
                                                        trio[1] = gap;
                                                        trio[2] = gap;

                                                        temp_state_111 = log(prev_state_to_100["111"]);
                                                        temp_state_110 = log(prev_state_to_100["110"]);
                                                        temp_state_101 = log(prev_state_to_100["101"]);
                                                        temp_state_100 = log(prev_state_to_100["100"]);
                                                        temp_state_011 = log(prev_state_to_100["11"]);
                                                        temp_state_010 = log(prev_state_to_100["10"]);
                                                        temp_state_001 = log(prev_state_to_100["1"]);

                                                        temp_emission1 = log(emission_prob[trio]);
                                                        temp_emission1 = log(emission_double[alpha_indexes_emissions[trio]][3]);

                                                        v100_1 = log(prev_state_to_111["100"]) + V111_b[dim1 + 1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                        v100_2 = log(prev_state_to_110["100"]) + V110_b[dim1 + 1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                        v100_3 = log(prev_state_to_101["100"]) + V101_b[dim1 + 1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                        v100_4 = log(prev_state_to_100["100"]) + V100_b[dim1 + 1][dim2][dim3] + log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                        v100_5 = log(prev_state_to_011["100"]) + V011_b[dim1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                        v100_6 = log(prev_state_to_010["100"]) + V010_b[dim1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                        v100_7 = log(prev_state_to_001["100"]) + V001_b[dim1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_001]][0]);

                                                        a = max_imum(v100_1, v100_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                                        b = v100_1 + v100_2 - a;
                                                        temp3 = a + log(1 + exp(b-a));
                                                        a = max_imum(v100_3, temp3, -DBL_MAX);
                                                        b = v100_3 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                                        a = max_imum(v100_4, temp3, -DBL_MAX);
                                                        b = v100_4 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                                        a = max_imum(v100_5, temp3, -DBL_MAX);
                                                        b = v100_5 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                                        a = max_imum(v100_6, temp3, -DBL_MAX);
                                                        b = v100_6 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                                        a = max_imum(v100_7, temp3, -DBL_MAX);
                                                        b = v100_7 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                                        temp_total = log(emission_prob[trio]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);

                                                        //V100[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);
                                                        //V100[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][3]) + maximum(v100_1, v100_2, v100_3, v100_4, v100_5, v100_6, v100_7);
                                                        //V100[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][3]) + temp3;
                                                        V100_b[dim1][dim2][dim3] = temp3;

                                                        // INS13 matrix

                                                        trio[0] = x1[dim1-1];
                                                        trio[1] = gap;
                                                        trio[2] = x3[dim3-1];

                                                        temp_state_111 = log(prev_state_to_101["111"]);
                                                        temp_state_110 = log(prev_state_to_101["110"]);
                                                        temp_state_101 = log(prev_state_to_101["101"]);
                                                        temp_state_100 = log(prev_state_to_101["100"]);
                                                        temp_state_011 = log(prev_state_to_101["11"]);
                                                        temp_state_010 = log(prev_state_to_101["10"]);
                                                        temp_state_001 = log(prev_state_to_101["1"]);

                                                        temp_emission1 = log(emission_prob[trio]);
                                                        temp_emission1 = log(emission_double[alpha_indexes_emissions[trio]][4]);

                                                        v101_1 = log(prev_state_to_111["101"]) + V111_b[dim1 + 1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                        v101_2 = log(prev_state_to_110["101"]) + V110_b[dim1 + 1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                        v101_3 = log(prev_state_to_101["101"]) + V101_b[dim1 + 1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                        v101_4 = log(prev_state_to_100["101"]) + V100_b[dim1 + 1][dim2][dim3] + log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                        v101_5 = log(prev_state_to_011["101"]) + V011_b[dim1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                        v101_6 = log(prev_state_to_010["101"]) + V010_b[dim1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                        v101_7 = log(prev_state_to_001["101"]) + V001_b[dim1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_001]][0]);

                                                        a = max_imum(v101_1, v101_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                                        b = v101_1 + v101_2 - a;
                                                        temp3 = a + log(1 + exp(b-a));
                                                        a = max_imum(v101_3, temp3, -DBL_MAX);
                                                        b = v101_3 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                                        a = max_imum(v101_4, temp3, -DBL_MAX);
                                                        b = v101_4 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                                        a = max_imum(v101_5, temp3, -DBL_MAX);
                                                        b = v101_5 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                                        a = max_imum(v101_6, temp3, -DBL_MAX);
                                                        b = v101_6 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                                        a = max_imum(v101_7, temp3, -DBL_MAX);
                                                        b = v101_7 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                                        temp_total = log(emission_prob[trio]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);

                                                        //V101[dim1][dim2][dim3] = log(emission_prob[trio]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);
                                                        //V101[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][4]) + maximum(v101_1, v101_2, v101_3, v101_4, v101_5, v101_6, v101_7);
                                                        //V101[dim1][dim2][dim3] = log(emission_double[alpha_indexes_emissions[trio]][4]) + temp3;
                                                        V101_b[dim1][dim2][dim3] = temp3;

                                                        // INS12 matrix

                                                        trio[0] = x1[dim1-1];
                                                        trio[1] = x2[dim2-1];
                                                        trio[2] = gap;

                                                        temp_state_111 = log(prev_state_to_110["111"]);
                                                        temp_state_110 = log(prev_state_to_110["110"]);
                                                        temp_state_101 = log(prev_state_to_110["101"]);
                                                        temp_state_100 = log(prev_state_to_110["100"]);
                                                        temp_state_011 = log(prev_state_to_110["11"]);
                                                        temp_state_010 = log(prev_state_to_110["10"]);
                                                        temp_state_001 = log(prev_state_to_110["1"]);

                                                        temp_emission1 = log(emission_prob[trio]);
                                                        temp_emission1 = log(emission_double[alpha_indexes_emissions[trio]][5]);

                                                        v110_1 = log(prev_state_to_111["110"]) + V111_b[dim1 + 1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_111]][6]);
                                                        v110_2 = log(prev_state_to_110["110"]) + V110_b[dim1 + 1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_110]][5]);
                                                        v110_3 = log(prev_state_to_101["110"]) + V101_b[dim1 + 1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_101]][4]);
                                                        v110_4 = log(prev_state_to_100["110"]) + V100_b[dim1 + 1][dim2][dim3] + log(emission_double[alpha_indexes_emissions[triplet_100]][3]);
                                                        v110_5 = log(prev_state_to_011["110"]) + V011_b[dim1][dim2 + 1][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_011]][2]);
                                                        v110_6 = log(prev_state_to_010["110"]) + V010_b[dim1][dim2 + 1][dim3] + log(emission_double[alpha_indexes_emissions[triplet_010]][1]);
                                                        v110_7 = log(prev_state_to_001["110"]) + V001_b[dim1][dim2][dim3 + 1] + log(emission_double[alpha_indexes_emissions[triplet_001]][0]);

                                                        a = max_imum(v110_1, v110_2, -DBL_MAX); // evaluating sum of logs for 'ALIGN' state.
                                                        b = v110_1 + v110_2 - a;
                                                        temp3 = a + log(1 + exp(b-a));
                                                        a = max_imum(v110_3, temp3, -DBL_MAX);
                                                        b = v110_3 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));                  // sum of v111_1, v111_2 and v111_3 = temp3
                                                        a = max_imum(v110_4, temp3, -DBL_MAX);
                                                        b = v110_4 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_4
                                                        a = max_imum(v110_5, temp3, -DBL_MAX);
                                                        b = v110_5 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_5
                                                        a = max_imum(v110_6, temp3, -DBL_MAX);
                                                        b = v110_6 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of temp3 and v111_6
                                                        a = max_imum(v110_7, temp3, -DBL_MAX);
                                                        b = v110_7 + temp3 - a;
                                                        temp3 = a + log(1 + exp(b-a));          // sum of all state transitions

                                                        temp_total = log(emission_prob[trio]) + maximum(v110_1, v110_2, v110_3, v110_4, v110_5, v110_6, v110_7);

                                                        V110_b[dim1][dim2][dim3] = temp3;
                                                }


                                        //}
                                }
                        }
                }

                /*// printing arrays
                cout << "V111_b: " << endl;
                for (dim1 = 0; dim1 <= x1.length()-1; dim1++)
                {
                        cout << "Dim1: " << dim1 << endl;
                        for (dim2 = 0; dim2 <= x2.length()-1; dim2++)
                        {
                                cout << "Dim2: " << dim2 << endl;
                                for (dim3 = 0; dim3 <= x3.length()-1; dim3++)
                                {
                                        cout << V111_b[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "V100_b: " << endl;
                for (dim1 = 0; dim1 <= x1.length()-1; dim1++)
                {
                        cout << "Dim1: " << dim1 << endl;
                        for (dim2 = 0; dim2 <= x2.length()-1; dim2++)
                        {
                                cout << "Dim2: " << dim2 << endl;
                                for (dim3 = 0; dim3 <= x3.length()-1; dim3++)
                                {
                                        cout << V100_b[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                cout << "V110_b: " << endl;
                for (dim1 = 0; dim1 <= x1.length()-1; dim1++)
                {
                        cout << "Dim1: " << dim1 << endl;
                        for (dim2 = 0; dim2 <= x2.length()-1; dim2++)
                        {
                                cout << "Dim2: " << dim2 << endl;
                                for (dim3 = 0; dim3 <= x3.length()-1; dim3++)
                                {
                                        cout << V110_b[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }*/
        }


        void max_expected_accuracy()
        {

                double trans_001, trans_010, trans_011, trans_100, trans_101, trans_110, trans_111, cost, temp, temp1, temp2;

                // BM array initializer
                for (int count1 = 0; count1 <= x1.length(); count1++)
                {
                        BM.push_back(rowvec1);
                        BM.push_back(rowvec1);
                        BM.push_back(rowvec1);
                        BM.push_back(rowvec1);
                        BM.push_back(rowvec1);
                        BM.push_back(rowvec1);
                        BM.push_back(rowvec1);

                        for (int count2 = 0; count2 <= x2.length(); count2++)
                        {
                                BM[count1].push_back(rowvec2);
                                BM[count1].push_back(rowvec2);
                                BM[count1].push_back(rowvec2);
                                BM[count1].push_back(rowvec2);
                                BM[count1].push_back(rowvec2);
                                BM[count1].push_back(rowvec2);
                                BM[count1].push_back(rowvec2);

                                for (int count3 = 0; count3 <= x3.length(); count3++)
                                {
                                        BM[count1][count2].push_back(0);
                                        BM[count1][count2].push_back(0);
                                        BM[count1][count2].push_back(0);
                                        BM[count1][count2].push_back(0);
                                        BM[count1][count2].push_back(0);
                                        BM[count1][count2].push_back(0);
                                        BM[count1][count2].push_back(0);
                                }
                        }
                }

                temp1 = exp(V111_f[1][1][1] + V111_b[0][0][0]);
                temp2 = exp(V111_f[2][2][2] + V111_b[1][1][1]);

                // matrix filling
                for (dim1 = 1; dim1 <= x1.length(); dim1++)
                {
                        for (dim2 = 1; dim2 <= x2.length(); dim2++)
                        {
                                for (dim3 = 1; dim3 <= x3.length(); dim3++)
                                {
                                        cost = exp(V111_f[dim1][dim2][dim3] + V111_b[dim1-1][dim2-1][dim3-1]);
                                        trans_111 = BM[dim1 - 1][dim2 - 1][dim3 - 1] + cost;
                                        trans_110 = BM[dim1 - 1][dim2 - 1][dim3];
                                        trans_101 = BM[dim1 - 1][dim2][dim3 - 1];
                                        trans_100 = BM[dim1 - 1][dim2][dim3];
                                        trans_011 = BM[dim1][dim2 - 1][dim3 - 1];
                                        trans_010 = BM[dim1][dim2 - 1][dim3];
                                        trans_001 = BM[dim1][dim2][dim3 - 1];
                                        temp = maximum(trans_001, trans_010, trans_011, trans_100, trans_101, trans_110, trans_111);

                                        BM[dim1][dim2][dim3] = maximum(trans_001, trans_010, trans_011, trans_100, trans_101, trans_110, trans_111);
                                }
                        }
                }

                // printing arrays
                for (dim1 = 1; dim1 <= x1.length(); dim1++)
                {
                        cout << "dim1: " << dim1 << endl;
                        for (dim2 = 1; dim2 <= x2.length(); dim2++)
                        {
                                cout << "dim2: " << dim2 << endl;
                                for (dim3 = 1; dim3 <= x3.length(); dim3++)
                                {
                                        //cout << BM[dim1][dim2][dim3] << " ";
                                }
                                cout << endl;
                        }
                }

                int wait = 1;

                // trace-back
                dim1 = x1.length();
                dim2 = x2.length();
                dim3 = x3.length();

                while (dim1 >= 1 && dim2 >= 1 && dim3 >= 1)
                {

                        cost = exp(V111_f[dim1][dim2][dim3] + V111_b[dim1-1][dim2-1][dim3-1]);

                        // aligned state
                        if (dim1 >= 1 && dim2 >= 1 && dim3 >= 1 && ((BM[dim1 - 1][dim2 - 1][dim3 - 1] + cost == BM[dim1][dim2][dim3])))
                        {
                                align1_reverse.push_back(x1[dim1-1]);
                                align2_reverse.push_back(x2[dim2-1]);
                                align3_reverse.push_back(x3[dim3-1]);
                                dim1--;
                                dim2--;
                                dim3--;
                        }

                        // insertion in seq 1,2 110
                        else if (dim1 >= 1 && dim2 >= 1 && ((BM[dim1 - 1][dim2 - 1][dim3] == BM[dim1][dim2][dim3])))
                        {
                                align1_reverse.push_back(x1[dim1-1]);
                                align2_reverse.push_back(x2[dim2-1]);
                                align3_reverse.push_back(gap);
                                dim1--;
                                dim2--;

                        }

                        // insertion in seq 1,3 101
                        else if (dim1 >= 1 && dim3 >= 1 && ((BM[dim1 - 1][dim2][dim3 - 1] == BM[dim1][dim2][dim3])))
                        {
                                align1_reverse.push_back(x1[dim1-1]);
                                align2_reverse.push_back(gap);
                                align3_reverse.push_back(x3[dim3-1]);
                                dim1--;
                                dim3--;
                        }

                        // insertion in seq 1, 100
                        else if (dim1 >= 1 && ((BM[dim1 - 1][dim2][dim3] == BM[dim1][dim2][dim3])))
                        {
                                align1_reverse.push_back(x1[dim1-1]);
                                align2_reverse.push_back(gap);
                                align3_reverse.push_back(gap);
                                dim1--;
                        }

                        // insertion in seq 2,3, 011
                        else if (dim2 >= 1 && dim3 >= 1 && ((BM[dim1][dim2 - 1][dim3 - 1] == BM[dim1][dim2][dim3])))
                        {
                                align1_reverse.push_back(gap);
                                align2_reverse.push_back(x2[dim2-1]);
                                align3_reverse.push_back(x3[dim3-1]);
                                dim2--;
                                dim3--;
                        }

                        // insertion in seq 2, 010
                        else if (dim2 >= 1 && ((BM[dim1][dim2 - 1][dim3] == BM[dim1][dim2][dim3])))
                        {
                                align1_reverse.push_back(gap);
                                align2_reverse.push_back(x2[dim2-1]);
                                align3_reverse.push_back(gap);
                                dim2--;
                        }

                        //      insertion in seq 3, 001
                        else //if (dim3 >= 1 && ((BM[dim1][dim2][dim3 - 1] == BM[dim1][dim2][dim3])))
                        {
                                align1_reverse.push_back(gap);
                                align2_reverse.push_back(gap);
                                align3_reverse.push_back(x3[dim3-1]);
                                //dim1--;
                                //dim2--;
                                dim3--;
                        }
                }

                // reversing arrays
                for (int i = align1_reverse.length(); i > 0; i--)
                {
                        align1.push_back(align1_reverse[i-1]);
                        align2.push_back(align2_reverse[i-1]);
                        align3.push_back(align3_reverse[i-1]);
                }


        }

};

#endif /* FUNCTIONS_H_ */