1 files changed, 47 insertions, 1 deletions

diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index a1c6c29..c82872e 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -494,7 +494,20 @@ void StandardizeVector(gsl_vector *y) {
 }
 
 // Calculate UtX (U gets transposed)
-// CalcUtX transforms a genotype matrix into the eigenspace by computing U^T × X, where U contains the eigenvectors from the kinship matrix decomposition. This transformation diagonalizes the correlation structure induced by relatedness, making subsequent likelihood calculations tractable and efficient. The function overwrites its input with the transformed result, using a temporary copy to avoid data corruption.
+/*
+CalcUtX transforms a genotype matrix by the eigenvector matrix from the kinship matrix decomposition.
+Purpose:
+
+Computes UtX = U^T × X, where:
+
+    U: Eigenvector matrix from the spectral decomposition of the kinship matrix K
+    X: Genotype matrix (individuals × SNPs)
+    UtX: Transformed genotype matrix
+
+This transformation rotates the genotype data into the eigenspace where the kinship-induced correlations become diagonal, simplifying likelihood calculations.
+
+CalcUtX transforms a genotype matrix into the eigenspace by computing U^T × X, where U contains the eigenvectors from the kinship matrix decomposition. This transformation diagonalizes the correlation structure induced by relatedness, making subsequent likelihood calculations tractable and efficient. The function overwrites its input with the transformed result, using a temporary copy to avoid data corruption.
+*/
 void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX) {
   gsl_matrix *X = gsl_matrix_safe_alloc(UtX->size1, UtX->size2);
   gsl_matrix_safe_memcpy(X, UtX);
@@ -632,3 +645,36 @@ double UcharToDouble02(const unsigned char c) { return (double)c * 0.01; }
 unsigned char Double02ToUchar(const double dosage) {
   return (int)(dosage * 100);
 }
+
+
+std::tuple<size_t, size_t, size_t> compute_ratio(size_t ni_total, const gsl_vector *gs) {
+  size_t n_0 = 0;
+  size_t n_1 = 0;
+  size_t n_2 = 0;
+  for (size_t i = 0; i < ni_total; ++i) { // read genotypes
+    double geno = gsl_vector_get(gs, i);
+    if (geno >= 0 && geno <= 0.5) {
+      n_0++;
+    }
+    if (geno > 0.5 && geno < 1.5) {
+      n_1++;
+    }
+    if (geno >= 1.5 && geno <= 2.0) {
+      n_2++;
+    }
+  }
+  return {n_0,n_1,n_2};
+}
+
+double compute_maf(size_t ni_total, size_t ni_test, size_t n_miss, const double *gs, const vector<int> &indicator) {
+  double maf = 0.0;
+
+  for (size_t i = 0; i < ni_total; ++i) { // read genotypes
+    if (indicator[i]) {
+      double geno = gs[i];
+      maf += geno; // 0..2 range
+    }
+  }
+  maf /= 2.0 * (double)(ni_test - n_miss); // Assumption is that geno value is between 0 and 2. FIXME
+  return maf;
+}