hmm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import time \n",
    "import json\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from pomegranate import State, DiscreteDistribution, HiddenMarkovModel\n",
    "from sklearn.model_selection import train_test_split\n",
    "from utils import *\n",
    "from Bio.Alphabet import IUPAC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GenerativeHMM(): \n",
    "    \n",
    "    def __init__(self, args, x_train=None, weights=None, verbose=True): \n",
    "        \"\"\"\n",
    "        Initializes the HMM to perform generative tasks\n",
    "        Parameters\n",
    "        ----------\n",
    "        args : dictionary\n",
    "            defines the hyper-parameters of the HMM\n",
    "        args.name : string \n",
    "            defines the name of the HMM\n",
    "        args.hidden_size : int \n",
    "            defines the hidden size\n",
    "        args.max_iterations: int\n",
    "            sets the max iterations\n",
    "        args.n_jobs: int\n",
    "            sets the number of cores to use\n",
    "        args.batch_size : int\n",
    "            sets the batch size\n",
    "        args.epochs : int \n",
    "            sets the epoch size \n",
    "        args.char_to_int : dict\n",
    "            a map from characters to index (integer) in the sequences\n",
    "        args.build_from_samples : boolean\n",
    "            build model from samples\n",
    "        \"\"\"\n",
    "        self.args = args\n",
    "        self.name = args[\"name\"]\n",
    "        self.hidden_size = args[\"hidden_size\"]\n",
    "        self.max_iterations = args[\"max_iterations\"]\n",
    "        self.n_jobs = args[\"n_jobs\"]\n",
    "        self.batch_size = args[\"batch_size\"]\n",
    "        self.epoch = args[\"epoch\"]\n",
    "        self.char_to_int = args[\"char_to_int\"]\n",
    "        self.vocabulary = [pair[0] for pair in sorted(self.char_to_int.items(), key = lambda x : x[1])]\n",
    "        self.indexes = [pair[1] for pair in sorted(self.char_to_int.items(), key = lambda x : x[1])]\n",
    "        self.emission_size = len(self.indexes)\n",
    "        if args[\"build_from_samples\"] and x_train is not None: \n",
    "            self.model = HiddenMarkovModel.from_samples(DiscreteDistribution, \n",
    "                                                    n_components = self.hidden_size, \n",
    "                                                    X = x_train, \n",
    "                                                    algorithm = 'baum-welch', \n",
    "                                                    return_history = True,\n",
    "                                                    verbose = verbose,\n",
    "                                                    max_iterations = self.max_iterations,\n",
    "                                                    n_jobs = self.n_jobs, \n",
    "                                                    weights = weights,\n",
    "                                                    #batch_size = self.batch_size,\n",
    "                                                    #batches_per_epoch = self.epochs\n",
    "                                               )[0]\n",
    "            \n",
    "        else: \n",
    "            self.build_model()\n",
    "        self.model.bake()\n",
    "\n",
    "    \n",
    "    def build_model(self): \n",
    "        distributions = []\n",
    "        for _ in range(self.hidden_size): \n",
    "            emission_probs = np.random.random(self.emission_size)\n",
    "            emission_probs = emission_probs / emission_probs.sum()\n",
    "            distributions.append(DiscreteDistribution(dict(zip(self.vocabulary, emission_probs))))\n",
    "        trans_mat = np.random.random((self.hidden_size, self.hidden_size))\n",
    "        trans_mat = trans_mat / trans_mat.sum(axis = 1, keepdims = 1)\n",
    "        starts = np.random.random((self.hidden_size))\n",
    "        starts = starts / starts.sum()\n",
    "        # testing initializations\n",
    "        np.testing.assert_almost_equal(starts.sum(), 1)\n",
    "        np.testing.assert_array_almost_equal(np.ones(self.hidden_size), trans_mat.sum(axis = 1))\n",
    "        self.model = HiddenMarkovModel.from_matrix(trans_mat, distributions, starts)\n",
    "        \n",
    "    def get_args(self): \n",
    "        return self.args\n",
    "        \n",
    "    def fit(self, x_train, weights=None, verbose=True):\n",
    "        \"\"\"\n",
    "        Fits the model on an HMM with self.hidden_size\n",
    "        \"\"\"    \n",
    "        return self.model.fit(x_train, \n",
    "                        algorithm = 'baum-welch', \n",
    "                        return_history = True, \n",
    "                        verbose = verbose,\n",
    "                        max_iterations = self.max_iterations,\n",
    "                        n_jobs = self.n_jobs, \n",
    "                        weights = weights,\n",
    "                        #batch_size = self.batch_size,\n",
    "                        #batches_per_epoch = self.epochs\n",
    "                   )\n",
    "    \n",
    "    def sample(self, n, length):\n",
    "        \"\"\"\n",
    "        Input:\n",
    "        n is number of samples\n",
    "        length is how long you want each sample to be\n",
    "        \"\"\"\n",
    "        return np.array([\"\".join(seq) for seq in self.model.sample(n = n, length = length)])\n",
    "            \n",
    "        \n",
    "    def predict(self, x_test): \n",
    "        \"\"\"\n",
    "        predict the log probability of obtaining the sequences in x_test\n",
    "        log(P(X1, X2, ..., X_test)) = sum(log(P(Xi)))\n",
    "        Input: x_test a list of sequences. should be 2 or 3 dimensional\n",
    "        \"\"\"\n",
    "        assert(len(np.array(x_test).shape) == 2 or len(np.array(x_test).shape) == 3)\n",
    "        return sum([self.model.log_probability(seq) for seq in np.array(x_test)])\n",
    "                \n",
    "    def show_model(self): \n",
    "        self.model.plot()\n",
    "        \n",
    "    def save_model(self, path): \n",
    "        with open(path, 'w') as f:\n",
    "            json.dump(self.model.to_json(), f)\n",
    "    \n",
    "    def load_model(self, path): \n",
    "        with open(path, 'r') as f:\n",
    "            json_model = json.load(f)\n",
    "        self.model = HiddenMarkovModel.from_json(json_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hmm_base_args(): \n",
    "    return {\n",
    "        \"name\" : \"base HMM\",\n",
    "        \"hidden_size\" : 5,\n",
    "        \"max_iterations\" : 10,\n",
    "        \"n_jobs\" : 1,\n",
    "        \"batch_size\" : 5,\n",
    "        \"epoch\" : 2,\n",
    "        \"char_to_int\" : {\"A\" : 0, \"C\" : 1, \"T\" : 2, \"G\" : 3},\n",
    "        \"build_from_samples\" : False\n",
    "    }\n",
    "\n",
    "def hmm_build_from_samples_args(): \n",
    "    args = hmm_base_args()\n",
    "    args[\"build_from_samples\"] = True\n",
    "    return args\n",
    "\n",
    "def hmm_amino_acid_args(): \n",
    "    args = hmm_base_args()\n",
    "    amino_acids = get_all_amino_acids()\n",
    "    indexes = list(range(len(amino_acids)))\n",
    "    assert(len(amino_acids) == 21)\n",
    "    assert(amino_acids == \"*\" + IUPAC.protein.letters) #*ACDEFGHIKLMNPQRSTVWY\n",
    "    args[\"char_to_int\"] = dict(zip(amino_acids, indexes))\n",
    "    return args\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_fit_dna_hmm(X_train_sequences): \n",
    "    args = hmm_base_args()\n",
    "    hmm = GenerativeHMM(args)\n",
    "    hmm.fit(X_train_sequences)\n",
    "    args[\"max_iterations\"] = 20\n",
    "    args[\"n_jobs\"] = 5\n",
    "    hmm = GenerativeHMM(args)\n",
    "    hmm.fit(X_train_sequences)\n",
    "    assert(hmm.char_to_int == {\"A\" : 0, \"C\" : 1, \"T\" : 2, \"G\" : 3})\n",
    "    assert(hmm.indexes == sorted(list(range(4))))\n",
    "    assert(hmm.vocabulary == list(\"ACTG\"))\n",
    "    \n",
    "def test_fit_amino_acid_hmm(X_train_sequences):\n",
    "    args = hmm_amino_acid_args()\n",
    "    hmm = GenerativeHMM(args)\n",
    "    hmm.fit(X_train_sequences)\n",
    "    args[\"max_iterations\"] = 20\n",
    "    args[\"n_jobs\"] = 5\n",
    "    hmm = GenerativeHMM(args)\n",
    "    hmm.fit(X_train_sequences)\n",
    "    amino_acid_alphabet = get_all_amino_acids()\n",
    "    assert(hmm.char_to_int == dict(zip(amino_acid_alphabet, list(range(len(amino_acid_alphabet))))))\n",
    "    assert(hmm.indexes == sorted(list(range(len(amino_acid_alphabet)))))\n",
    "    assert(hmm.vocabulary == list(amino_acid_alphabet))\n",
    "    \n",
    "def test_sample_and_predict_dna_hmm(x_train_sequences): \n",
    "    args = hmm_base_args()\n",
    "    hmm = GenerativeHMM(args)\n",
    "    hmm.fit(x_train_sequences, verbose=False)\n",
    "    seq1, seq2 = tuple(hmm.sample(2, 714))\n",
    "    np.testing.assert_almost_equal(hmm.model.probability(seq1), np.e ** hmm.predict([list(seq1)]))\n",
    "    total = 0\n",
    "    for i in \"ACTG\": \n",
    "        for j in \"ACTG\": \n",
    "            for k in \"ACTG\":\n",
    "                codon = i + j + k\n",
    "                np.testing.assert_almost_equal(hmm.model.probability(codon), np.e ** hmm.predict([list(codon)]))\n",
    "                total += np.e ** hmm.predict([list(codon)])\n",
    "    np.testing.assert_almost_equal(1, total)\n",
    "    \n",
    "def test_sample_and_predict_amino_acid_hmm(x_train_sequences): \n",
    "    args = hmm_amino_acid_args()\n",
    "    hmm = GenerativeHMM(args)\n",
    "    hmm.fit(x_train_sequences, verbose=False)\n",
    "    wild_type_amino_acid = get_wild_type_amino_acid_sequence()\n",
    "    seq1, seq2 = tuple(hmm.sample(2, len(wild_type_amino_acid)))\n",
    "    print(\"{0} amino acids away from wild type!\".format(count_substring_mismatch(seq1, \n",
    "                                                            wild_type_amino_acid)))\n",
    "    np.testing.assert_almost_equal(hmm.model.probability(seq1), np.e ** hmm.predict([list(seq1)]))\n",
    "    total = 0\n",
    "    amino_acid_alphabet = get_all_amino_acids()\n",
    "    for i in amino_acid_alphabet: \n",
    "        for j in amino_acid_alphabet: \n",
    "            amino_acid = i + j\n",
    "            np.testing.assert_almost_equal(hmm.model.probability(amino_acid), np.e ** hmm.predict([list(amino_acid)]))\n",
    "            total += np.e ** hmm.predict([list(amino_acid)])\n",
    "    np.testing.assert_almost_equal(total, 1)   \n",
    "    \n",
    "def test_fit_dna_hmm_from_samples(X_train_sequences): \n",
    "    args = hmm_build_from_samples_args()\n",
    "    hmm = GenerativeHMM(args, X_train_sequences)\n",
    "    #test max iterations and n_jobs \n",
    "    args[\"max_iterations\"] = 20\n",
    "    args[\"n_jobs\"] = 5\n",
    "    hmm = GenerativeHMM(args, X_train_sequences)\n",
    "    assert(hmm.char_to_int == {\"A\" : 0, \"C\" : 1, \"T\" : 2, \"G\" : 3})\n",
    "    assert(hmm.indexes == sorted(list(range(4))))\n",
    "    assert(hmm.vocabulary == list(\"ACTG\"))\n",
    "\n",
    "def test_fit_dna_hmm_weights(X_train_sequences):\n",
    "    args = hmm_base_args()\n",
    "    weights = np.identity(4)\n",
    "    weights = np.vstack([weights, [0.25, 0.25, 0.25, 0.25]])\n",
    "    for weight in weights:\n",
    "        counts = {\"A\" : 0, \"C\" : 0, \"T\" : 0, \"G\" : 0}\n",
    "        hmm = GenerativeHMM(args)\n",
    "        hmm.fit(X_train_sequences, weight, verbose=False)\n",
    "        json_model = json.loads(hmm.model.to_json())\n",
    "        for state in json_model[\"states\"]: \n",
    "            if state is not None and state[\"distribution\"] is not None:\n",
    "                mp = state[\"distribution\"][\"parameters\"][0]\n",
    "                for k, v in mp.items(): \n",
    "                    counts[k] = counts[k] + v\n",
    "        print(\"Weights:\", weight, \"\\nCounts:\", counts)    \n",
    "        \n",
    "\n",
    "def test_fit_dna_hmm_from_samples_weights(X_train_sequences):\n",
    "    args = hmm_build_from_samples_args()\n",
    "    weights = np.identity(4)\n",
    "    weights = np.vstack([weights, [0.25, 0.25, 0.25, 0.25]])\n",
    "    for weight in weights: \n",
    "        counts = {\"A\" : 0, \"C\" : 0, \"T\" : 0, \"G\" : 0}\n",
    "        hmm = GenerativeHMM(args, X_train_sequences, weight, verbose=False)\n",
    "        json_model = json.loads(hmm.model.to_json())\n",
    "        for state in json_model[\"states\"]: \n",
    "            if state is not None and state[\"distribution\"] is not None:\n",
    "                mp = state[\"distribution\"][\"parameters\"][0]\n",
    "                for k, v in mp.items(): \n",
    "                    counts[k] = counts[k] + v\n",
    "        print(\"Weights:\", weight, \"\\nCounts:\", counts) \n",
    "        \n",
    "def test_save_and_load_hmm(X_train_sequences, amino_acid = True): \n",
    "    if amino_acid:\n",
    "        args = hmm_amino_acid_args()\n",
    "    else:\n",
    "        args = hmm_base_args()\n",
    "    hmm = GenerativeHMM(args)\n",
    "    hmm.fit(X_train_sequences, verbose=False)\n",
    "    hmm.save_model(\"./models/test.json\")\n",
    "    cached_hmm = GenerativeHMM(args)\n",
    "    cached_hmm.load_model(\"./models/test.json\")\n",
    "    for i in \"ACTG\": \n",
    "        for j in \"ACTG\": \n",
    "            for k in \"ACTG\":\n",
    "                codon = i + j + k\n",
    "                np.testing.assert_almost_equal(hmm.predict([list(codon)]), cached_hmm.predict([list(codon)]))\n",
    "\n",
    "def test_args(X_train_sequences, amino_acid=True):\n",
    "    if amino_acid:\n",
    "        args = hmm_amino_acid_args()\n",
    "    else:\n",
    "        args = hmm_base_args()\n",
    "    hmm = GenerativeHMM(args)\n",
    "    hmm.fit(X_train_sequences, verbose=False)\n",
    "    assert(hmm.args == args and hmm.get_args() == args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Improvement: 4809.994809303782\tTime (s): 0.07913\n",
      "[2] Improvement: 195.9385597890505\tTime (s): 0.0757\n",
      "[3] Improvement: 208.48861080520146\tTime (s): 0.07947\n",
      "[4] Improvement: 234.98140145662182\tTime (s): 0.09808\n",
      "[5] Improvement: 267.5174148558872\tTime (s): 0.08979\n",
      "[6] Improvement: 297.23517756792717\tTime (s): 0.08581\n",
      "[7] Improvement: 313.36477050112444\tTime (s): 0.07807\n",
      "[8] Improvement: 311.474847865029\tTime (s): 0.07326\n",
      "[9] Improvement: 298.2946238792065\tTime (s): 0.07482\n",
      "[10] Improvement: 282.7325340189418\tTime (s): 0.07502\n",
      "Total Training Improvement: 7220.022750042772\n",
      "Total Training Time (s): 0.9153\n",
      "[1] Improvement: 5496.2382646826445\tTime (s): 0.03234\n",
      "[2] Improvement: 171.56578296980297\tTime (s): 0.03409\n",
      "[3] Improvement: 205.562677476948\tTime (s): 0.03475\n",
      "[4] Improvement: 261.3211669871089\tTime (s): 0.03227\n",
      "[5] Improvement: 333.1231366408465\tTime (s): 0.03396\n",
      "[6] Improvement: 403.11317243370286\tTime (s): 0.04896\n",
      "[7] Improvement: 439.38394473800145\tTime (s): 0.0462\n",
      "[8] Improvement: 431.51644405836123\tTime (s): 0.06008\n",
      "[9] Improvement: 399.83584396974766\tTime (s): 0.04492\n",
      "[10] Improvement: 363.1435449078708\tTime (s): 0.06252\n",
      "[11] Improvement: 330.6027158476354\tTime (s): 0.03629\n",
      "[12] Improvement: 327.6692938307242\tTime (s): 0.03662\n",
      "[13] Improvement: 389.228791632042\tTime (s): 0.03741\n",
      "[14] Improvement: 491.5858300949694\tTime (s): 0.03409\n",
      "[15] Improvement: 522.9635665175592\tTime (s): 0.0326\n",
      "[16] Improvement: 425.02317056605534\tTime (s): 0.03494\n",
      "[17] Improvement: 274.13380143568793\tTime (s): 0.05058\n",
      "[18] Improvement: 174.75386209780117\tTime (s): 0.05032\n",
      "[19] Improvement: 143.84191594095319\tTime (s): 0.05225\n",
      "[20] Improvement: 158.95840472975397\tTime (s): 0.05503\n",
      "Total Training Improvement: 11743.565331558217\n",
      "Total Training Time (s): 0.9493\n",
      "228 amino acids away from wild type!\n"
     ]
    }
   ],
   "source": [
    "test_args(test_amino_acid_sequences)\n",
    "test_fit_amino_acid_hmm(test_amino_acid_sequences)\n",
    "test_sample_and_predict_amino_acid_hmm(test_amino_acid_sequences)\n",
    "test_save_and_load_hmm(test_amino_acid_sequences)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Improvement: 3936.230649036137\tTime (s): 0.2099\n",
      "[2] Improvement: 149.12637396418722\tTime (s): 0.2391\n",
      "[3] Improvement: 108.95387729156937\tTime (s): 0.2187\n",
      "[4] Improvement: 85.98553325257672\tTime (s): 0.2246\n",
      "[5] Improvement: 72.22916251717834\tTime (s): 0.2097\n",
      "[6] Improvement: 63.96028205938637\tTime (s): 0.2168\n",
      "[7] Improvement: 59.253396770436666\tTime (s): 0.2162\n",
      "[8] Improvement: 57.03780027135508\tTime (s): 0.2401\n",
      "[9] Improvement: 56.68740225168585\tTime (s): 0.2234\n",
      "[10] Improvement: 57.828405017979094\tTime (s): 0.2073\n",
      "Total Training Improvement: 4647.292882432492\n",
      "Total Training Time (s): 2.4972\n",
      "[1] Improvement: 3721.3673773265473\tTime (s): 0.07131\n",
      "[2] Improvement: 22.49783210310852\tTime (s): 0.0714\n",
      "[3] Improvement: 20.538204476280953\tTime (s): 0.07369\n",
      "[4] Improvement: 19.737539391193422\tTime (s): 0.09233\n",
      "[5] Improvement: 19.73047075814975\tTime (s): 0.07716\n",
      "[6] Improvement: 20.38657226327632\tTime (s): 0.07107\n",
      "[7] Improvement: 21.653334047354292\tTime (s): 0.08918\n",
      "[8] Improvement: 23.525946802008548\tTime (s): 0.1207\n",
      "[9] Improvement: 26.044249900471186\tTime (s): 0.1051\n",
      "[10] Improvement: 29.295794634876074\tTime (s): 0.07824\n",
      "[11] Improvement: 33.42400787392398\tTime (s): 0.08292\n",
      "[12] Improvement: 38.64407945089624\tTime (s): 0.08803\n",
      "[13] Improvement: 45.27048245223705\tTime (s): 0.08207\n",
      "[14] Improvement: 53.76226901212067\tTime (s): 0.07948\n",
      "[15] Improvement: 64.79652235706453\tTime (s): 0.1041\n",
      "[16] Improvement: 79.38794142623374\tTime (s): 0.1065\n",
      "[17] Improvement: 99.0860321851651\tTime (s): 0.1169\n",
      "[18] Improvement: 126.30517463706201\tTime (s): 0.08921\n",
      "[19] Improvement: 164.88293333542242\tTime (s): 0.07933\n",
      "[20] Improvement: 221.0194192546478\tTime (s): 0.09964\n",
      "Total Training Improvement: 4851.35618368804\n",
      "Total Training Time (s): 2.0156\n",
      "[1] Improvement: 4923.43666946009\tTime (s): 0.2148\n",
      "[2] Improvement: 23.526267938039382\tTime (s): 0.2144\n",
      "[3] Improvement: 20.35255074854649\tTime (s): 0.2137\n",
      "[4] Improvement: 17.053787234297488\tTime (s): 0.2054\n",
      "[5] Improvement: 14.551319871359738\tTime (s): 0.2293\n",
      "[6] Improvement: 13.093626063142437\tTime (s): 0.2077\n",
      "[7] Improvement: 12.437872376816813\tTime (s): 0.2129\n",
      "[8] Improvement: 12.2931049913459\tTime (s): 0.2045\n",
      "[9] Improvement: 12.458137739347876\tTime (s): 0.1999\n",
      "[10] Improvement: 12.817085512899212\tTime (s): 0.1948\n",
      "Total Training Improvement: 5062.020421935886\n",
      "Total Training Time (s): 2.3856\n",
      "[1] Improvement: 10409.974112110387\tTime (s): 0.0739\n",
      "[2] Improvement: 54.757567181717604\tTime (s): 0.08188\n",
      "[3] Improvement: 40.912123348229215\tTime (s): 0.07485\n",
      "[4] Improvement: 31.34543094894616\tTime (s): 0.1164\n",
      "[5] Improvement: 25.504439059266588\tTime (s): 0.0957\n",
      "[6] Improvement: 21.827727430936648\tTime (s): 0.07923\n",
      "[7] Improvement: 19.31719153409358\tTime (s): 0.1153\n",
      "[8] Improvement: 17.480315767999855\tTime (s): 0.07884\n",
      "[9] Improvement: 16.083993895706953\tTime (s): 0.08466\n",
      "[10] Improvement: 15.013837267702911\tTime (s): 0.1265\n",
      "[11] Improvement: 14.210333670765976\tTime (s): 0.1236\n",
      "[12] Improvement: 13.640709023995441\tTime (s): 0.09262\n",
      "[13] Improvement: 13.286062233892153\tTime (s): 0.1583\n",
      "[14] Improvement: 13.135320051689632\tTime (s): 0.1268\n",
      "[15] Improvement: 13.182545963354642\tTime (s): 0.1008\n",
      "[16] Improvement: 13.42602546676062\tTime (s): 0.09778\n",
      "[17] Improvement: 13.868284571522963\tTime (s): 0.1011\n",
      "[18] Improvement: 14.516535364033189\tTime (s): 0.0815\n",
      "[19] Improvement: 15.383237222747994\tTime (s): 0.0754\n",
      "[20] Improvement: 16.486588657964603\tTime (s): 0.1044\n",
      "Total Training Improvement: 10793.352380771714\n",
      "Total Training Time (s): 2.2121\n",
      "Weights: [1. 0. 0. 0.] \n",
      "Counts: {'A': 3.2366628269343174, 'C': 0.0, 'T': 1.7633371730656822, 'G': 0.0}\n",
      "Weights: [0. 1. 0. 0.] \n",
      "Counts: {'A': 0.0, 'C': 3.6432603894713087, 'T': 0.0, 'G': 1.3567396105286909}\n",
      "Weights: [0. 0. 1. 0.] \n",
      "Counts: {'A': 1.6146422698473648, 'C': 0.0, 'T': 3.385357730152635, 'G': 0.0}\n",
      "Weights: [0. 0. 0. 1.] \n",
      "Counts: {'A': 0.0, 'C': 1.3789836122310095, 'T': 0.0, 'G': 3.6210163877689907}\n",
      "Weights: [0.25 0.25 0.25 0.25] \n",
      "Counts: {'A': 1.051738312366067, 'C': 1.1586217545660378, 'T': 1.2397422279527461, 'G': 1.5498977051151492}\n",
      "Weights: [1. 0. 0. 0.] \n",
      "Counts: {'A': 3.6639529160844173, 'C': 0.0, 'T': 1.3360470839155827, 'G': 0.0}\n",
      "Weights: [0. 1. 0. 0.] \n",
      "Counts: {'A': 0.0, 'C': 3.4747403813734152, 'T': 0.0, 'G': 1.5252596186265845}\n",
      "Weights: [0. 0. 1. 0.] \n",
      "Counts: {'A': 1.613648565579176, 'C': 0.0, 'T': 3.386351434420824, 'G': 0.0}\n",
      "Weights: [0. 0. 0. 1.] \n",
      "Counts: {'A': 0.0, 'C': 1.607338867059708, 'T': 0.0, 'G': 3.392661132940292}\n",
      "Weights: [0.25 0.25 0.25 0.25] \n",
      "Counts: {'A': 2.1551585661785237, 'C': 1.1984595425881228, 'T': 0.8447389753696606, 'G': 0.8016429158636922}\n"
     ]
    }
   ],
   "source": [
    "test_fit_dna_hmm(test_dna_sequences)\n",
    "test_sample_and_predict_dna_hmm(test_dna_sequences)\n",
    "test_fit_dna_hmm_from_samples(test_dna_sequences)\n",
    "test_save_and_load_hmm(test_dna_sequences)\n",
    "\n",
    "synthetic_data = np.array([[\"A\", \"A\", \"A\", \"T\"], [\"C\", \"C\", \"C\", \"G\"], [\"T\", \"T\", \"T\", \"A\"], [\"G\", \"G\", \"G\", \"C\"]])\n",
    "test_fit_dna_hmm_weights(synthetic_data)\n",
    "test_fit_dna_hmm_from_samples_weights(synthetic_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading data...\n",
      "Finished loading data in 4.46 seconds\n"
     ]
    }
   ],
   "source": [
    "print(\"Loading data...\")\n",
    "start_time = time.time()\n",
    "X_train, X_test, y_train, y_test = load_gfp_data(\"./data/gfp_dna_\") #./data/gfp_dna_\n",
    "mutated_df = load_saved_mutated_gfp_data()\n",
    "assert(X_train[0] == get_wild_type_dna_sequence())\n",
    "assert(count_substring_mismatch(X_train[999], get_wild_type_dna_sequence()) == 3)\n",
    "print(\"Finished loading data in {0:.2f} seconds\".format(time.time() - start_time))\n",
    "\n",
    "base_args = hmm_base_args()\n",
    "build_from_samples_args = hmm_build_from_samples_args()\n",
    "test_dna_sequences = np.array([list(seq) for seq in X_train[0:100]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading data...\n",
      "Finished loading data in 1.64 seconds\n"
     ]
    }
   ],
   "source": [
    "print(\"Loading data...\")\n",
    "start_time = time.time()\n",
    "X_train, X_test, y_train, y_test = load_gfp_data(\"./data/gfp_amino_acid_\") #./data/gfp_amino_acid_\n",
    "assert(X_train[0] == get_wild_type_amino_acid_sequence())\n",
    "assert(count_substring_mismatch(X_train[999], get_wild_type_amino_acid_sequence()) == 3)\n",
    "print(\"Finished loading data in {0:.2f} seconds\".format(time.time() - start_time))\n",
    "\n",
    "base_amino_acid_args = hmm_amino_acid_args()\n",
    "test_amino_acid_sequences = np.array([list(seq) for seq in X_train[0:100]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}