 }
 
 
-void getFreqSingle(int* alignment, float* f_single, size_t B, size_t N){
+void getFreqSingle(int* aln, float* f_single, size_t B, size_t N){
 
     for (int i = 0; i < N; i++){
+        /*
         for (int k = 0; k < Q; k++){
             float f_ki = 0.0;
             for (int b = 0; b < B; b++){
-                f_ki += 1/(float)B * delta(alignment[b*N + i], k+1);
+                f_ki += 1/(float)B * delta(aln[b*N + i], k+1);
             }
             f_single[k*N + i] = f_ki;
-        }
-    }
-}
-
-__global__ void getFreqSingleOnDevice(int* alignment, float* f_single, size_t B, size_t N){
+        }*/
+        for (int k = 0; k < Q; k++)
+            f_single[k*N + i] = 0.0;
+        
+        for (int b = 0; b < B; b++){
+            int k = aln[b*N + i] - 1;
+            f_single[k*N + i] += 1.0; 
+        }
+
+        for (int k = 0; k < Q; k++)
+            f_single[k*N + i] /= (float)B;
+    }
+}
+
+__global__ void getFreqSingleOnDevice(int* aln, float* f_single, size_t B, size_t N, int Q){
 
     // Grid dimensions: N x Q
     int i = threadIdx.x;
-    int k = threadIdx.y;
-
+    //int k = threadIdx.y;
+
+    /*
     float f_ki = 0.0;
     for (int b = 0; b < B; b++){
-        if (alignment[b*N + i] == k+1){
+        if (aln[b*N + i] == k+1){
             f_ki += 1/(float)B;
         }
     }
     f_single[k*N + i] = f_ki;
-}
-
-
-void getFreqPair(int* alignment, float* f_pair, size_t B, size_t N){
-
+    */
+    if(i < N){
+    for (int k = 0; k < Q; k++)
+        f_single[k*N + i] = 0.0;
+    
+    for (int b = 0; b < B; b++){
+        int k = aln[b*N + i] - 1;
+        f_single[k*N + i] += 1.0; 
+    }
+
+    for (int k = 0; k < Q; k++)
+        f_single[k*N + i] /= (float)B;
+
+    }
+}
+
+
+void getFreqPair(int* aln, float* f_pair, size_t B, size_t N){
+
+    /*
     for (int i = 0; i < N; i++){
     for (int j = 0; j < N; j++){
         for (int k = 0; k < Q; k++){
         for (int l = 0; l < Q; l++){
             float f_klij = 0.0;
             for (int b = 0; b < B; b++){
-                f_klij += 1/(float)B * delta(alignment[b*N + i], k+1) * delta(alignment[b*N + j], l+1);
+                f_klij += 1/(float)B * delta(aln[b*N + i], k+1) * delta(aln[b*N + j], l+1);
             }
             f_pair[(Q*i + k) * Q*N + (Q*j + l)] = f_klij;
         }
         }
     }
     }
-}
-
-__global__ void getFreqPairOnDevice(int* alignment, float* f_pair, size_t B, size_t N, int Q){
+    */
+
+    for (int i = 0; i < N; i++){
+    for (int j = 0; j < N; j++){
+        for (int k = 0; k < Q; k++)
+        for (int l = 0; l < Q; l++)
+            f_pair[(Q*i + k) * Q*N + (Q*j + l)] = 0.0;
+
+        for (int b = 0; b < B; b++){
+            int k = aln[b*N + i] - 1;
+            int l = aln[b*N + j] - 1;
+            f_pair[(Q*i + k) * Q*N + (Q*j + l)] += 1.0;
+        }
+
+        for (int k = 0; k < Q; k++){
+        for (int l = 0; l < Q; l++){
+            f_pair[(Q*i + k) * Q*N + (Q*j + l)] /= (float)B;
+            //f_pair[(Q*j + l) * Q*N + (Q*i + k)] = f_pair[(Q*i + k) * Q*N + (Q*j + l)];
+        }
+        }
+    }
+    }
+
+}
+
+__global__ void getFreqPairOnDevice(int* aln, float* f_pair, size_t B, size_t N, int Q){
 
     // Grid dimensions: N x N
     int i = threadIdx.x;
     int j = threadIdx.y;
 
+    /*
     for (int k = 0; k < Q; k++){
     for (int l = 0; l < Q; l++){
         float f_klij = 0.0;
         for (int b = 0; b < B; b++){
-            if ((alignment[b*N + i] == k+1) && (alignment[b*N + j] == l+1)){
+            if ((aln[b*N + i] == k+1) && (aln[b*N + j] == l+1)){
                 f_klij += 1/(float)B;
             }
         }
         f_pair[(Q*i + k) * Q*N + (Q*j + l)] = f_klij;
     }
     }
-}
-
+    */
+    if(i<N && j<N){
+    for (int k = 0; k < Q; k++)
+    for (int l = 0; l < Q; l++)
+        f_pair[(Q*i + k) * Q*N + (Q*j + l)] = 0.0;
+
+    for (int b = 0; b < B; b++){
+        int k = aln[b*N + i] - 1;
+        int l = aln[b*N + j] - 1;
+        f_pair[(Q*i + k) * Q*N + (Q*j + l)] += 1.0;
+    }
+
+    for (int k = 0; k < Q; k++){
+    for (int l = 0; l < Q; l++){
+        f_pair[(Q*i + k) * Q*N + (Q*j + l)] /= (float)B;
+        //f_pair[(Q*j + l) * Q*N + (Q*i + k)] = f_pair[(Q*i + k) * Q*N + (Q*j + l)];
+    }
+    }
+    }
+}
+    
 
 void getCovMat(float* f_single, float* f_pair, float* cov_mat, size_t B, size_t N){
 


    // calculate column-wise amino acid frequencies
    float* f_single = (float*) std::malloc(Q*N * sizeof(float));
    getFreqSingle(alignment, f_single, B, N);

    // calculate column-wise amino acid frequencies
    // for each possible pair of amino acids and columns
    float* f_pair = (float*) std::malloc(Q*N * Q*N * sizeof(float));
    getFreqPair(alignment, f_pair, B, N);

    // calculate covariance matrix from frequencies
    float* cov_mat = (float*) std::malloc(Q*N * Q*N * sizeof(float));

#include <map>
#include <cstdlib>

#include <stdio.h> 
#include <stdlib.h> 
#include <stdarg.h> 
#include <string.h> 
#include <ctype.h> 
#include <math.h> 
#include <unistd.h> 
 
#include <time.h>
#include <assert.h>

