TRIPLE ALIGNMENT!!

nifets · May 18, 2023 · 4cc76d8 · 4cc76d8
1 parent 9aa5142
commit 4cc76d8
Show file tree

Hide file tree

Showing 10 changed files with 241 additions and 1,281 deletions.
diff --git a/.ipynb_checkpoints/Ancestor Reconstruction-checkpoint.ipynb b/.ipynb_checkpoints/Ancestor Reconstruction-checkpoint.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 4,
    "id": "c79e74a3",
    "metadata": {},
    "outputs": [
@@ -17,10 +17,36 @@
     {
      "data": {
       "text/plain": [
-       "-1281.439237586948"
+       "142×146 Matrix{Float64}:\n",
+       " -25.0358  -24.7162  -25.0068  -33.7361   …  -28.0232  -66.9887   -13.434\n",
+       " -36.1451  -15.4942  -16.2317  -25.3725      -19.0279  -32.9741   -13.6597\n",
+       " -33.6691  -21.2741  -14.5955  -25.4042      -20.5477  -44.5028   -11.7251\n",
+       " -39.7601  -34.1013  -33.2484   -5.02986     -35.2034  -30.3247    -3.88709\n",
+       " -40.4012  -32.9982  -32.0806   -7.61196     -34.7408  -27.4761    -3.17999\n",
+       " -39.8164  -35.0303  -33.3116   -9.24855  …  -34.5567  -28.6569    -3.61262\n",
+       " -39.9721  -33.9442  -32.4771   -8.79749     -35.2304  -27.4783    -3.5067\n",
+       " -40.3626  -33.7294  -29.2066   -8.62386     -34.5341  -29.1127    -3.54312\n",
+       " -39.4897  -34.8798  -32.3605   -9.88389     -33.899   -27.2323    -3.96756\n",
+       " -42.5864  -31.472   -31.3898   -9.33567     -33.4383  -29.8048    -3.36862\n",
+       " -39.4256  -33.5821  -32.9801   -8.75491  …  -35.4168  -27.1953    -3.52709\n",
+       " -41.3326  -33.0187  -31.1069   -7.76237     -33.8885  -28.4409    -3.16163\n",
+       " -41.5296  -33.063   -31.1095   -7.71654     -33.9511  -28.0034    -3.15441\n",
+       "   ⋮                                      ⋱                         ⋮\n",
+       " -40.9774  -33.633   -31.1651   -7.74883  …  -33.7832  -26.5294    -3.37743\n",
+       " -41.7462  -31.8062  -31.7774   -9.38849     -33.8059  -30.4752    -3.39513\n",
+       " -39.0229  -33.8772  -31.5923   -8.10748     -33.8955  -29.7599    -3.47554\n",
+       " -41.0386  -33.4216  -28.7636   -8.64531     -34.0483  -29.0397    -3.52754\n",
+       " -41.1785  -32.1647  -32.3694   -9.24471     -34.5342  -29.6882    -3.38969\n",
+       " -37.5481  -26.5371  -35.3924   -8.84838  …  -32.1875  -27.3223    -3.38098\n",
+       " -43.4024  -27.7105  -26.5937  -10.1882      -29.1568  -19.1809    -4.9569\n",
+       " -37.2676  -30.633   -33.4871   -7.85298     -31.7102  -26.8061    -3.61648\n",
+       " -41.8893  -26.5085  -27.5764  -11.1675      -28.0125  -18.6169    -5.37511\n",
+       " -37.8765  -32.0577  -35.9568  -10.1224      -29.0982  -25.5581    -4.17535\n",
+       " -60.3151  -35.7187  -37.3722  -17.2967   …  -37.0523  -13.4074    -8.95197\n",
+       " -18.4819  -14.2604  -10.9914   -3.84752     -14.249    -9.21093    0.0"
       ]
      },
-     "execution_count": 19,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -54,8 +80,8 @@
     "λ = 0.015; μ = 0.0156; r = 0.6\n",
     "align_model = TKF92([t], λ, μ, r; known_ancestor=true)\n",
     "pair_hmm = PairDataHMM(align_model, num_sites(Y), num_sites(Z));\n",
-    "emission_lps = fulljointlogpdf(ξ, t, Y, Z);\n",
-    "logpdf(pair_hmm, emission_lps)"
+    "pair_chain_dist = ChainJointDistribution(ξ, align_model)\n",
+    "emission_lps = fulljointlogpdf(ξ, t, Y, Z)"
    ]
   },
   {
@@ -76,6 +102,7 @@
     }
    ],
    "source": [
+    "# compute emission lp for each site in alignment, given pairwise emission matrix\n",
     "function logpdfs(alignment, emission_lps)\n",
     "    res = Real[]\n",
     "    N, M = size(emission_lps) .- 1\n",
@@ -97,6 +124,7 @@
     "    return res\n",
     "end\n",
     "\n",
+    "# compute emission lps for specific regime choices per site\n",
     "function logpdfsregime(alignment, regimes, Y, Z, ξ, t)\n",
     "    M = length(alignment)\n",
     "    C = num_coords(ξ)\n",
@@ -319,18 +347,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "function ancestor_sampling(M_YZ, Y, Z, t, ξ)\n",
+    "function ancestor_sampling(M_YZ::Alignment, Y::ObservedChain, Z:;ObservedChain, \n",
+    "                           t::Real, ξ::MixtureProductProcess)\n",
     "    # step 1 - sample XYZ alignment\n",
     "    M_XYZ = sample_anc_alignment(M_YZ, Y, Z, t, ξ)\n",
     "    \n",
     "    # step 2 - sample coordinates of X\n",
     "    alignment = M_XYZ\n",
     "    X_mask = mask(alignment, [[1], [0,1], [0,1]])\n",
-    "    alignmentX = alignment[X_mask]\n",
+    "    alignmentX = slice(alignment, X_mask)\n",
     "    Y_mask = mask(alignment, [[0,1], [1], [0,1]])\n",
-    "    alignmentY = alignment[Y_mask]\n",
+    "    alignmentY = slice(alignment, Y_mask)\n",
     "    Z_mask = mask(alignment, [[0,1], [0,1], [1]])\n",
-    "    alignmentZ = alignment[Z_mask]\n",
+    "    alignmentZ = slice(alignment, Z_mask)\n",
     "    \n",
     "    M = length(alignment)\n",
     "    regimes = ones(M)\n",
@@ -386,7 +415,7 @@
     "        dataX111 .= sample_anc_coords(p, [dataY111, dataZ111], t)\n",
     "    end\n",
     "    \n",
-    "    X = ObservedData(dataX)\n",
+    "    X = ObservedChain(dataX)\n",
     "end\n",
     "\n",
     "function trajectory_reconstruction(Y, Z, align_model, ξ, t;levels=1)\n",