Untitled.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "As= [[-1.53240875e-07  6.07056185e-08  3.61990772e-08  5.63361792e-08]\n",
      " [ 7.26517041e-08 -8.75254095e-07  5.80768946e-08  7.44525496e-07]\n",
      " [ 4.64970906e-08  6.23325291e-08 -2.15517146e+00  2.15517136e+00]\n",
      " [ 7.59646218e-08  8.38854519e-07  2.26244233e+00 -2.26244324e+00]]\n",
      "P = [[9.99999998e-01 6.07056184e-10 3.64426875e-10 5.60925689e-10]\n",
      " [7.26517040e-10 9.99999991e-01 6.57597477e-10 7.36842636e-09]\n",
      " [4.68100038e-10 7.05783500e-10 9.78917388e-01 2.10826108e-02]\n",
      " [7.56361336e-10 8.30198265e-09 2.21319715e-02 9.77868019e-01]]\n",
      "[[9.99999951e-01 1.94257971e-08 1.32479474e-08 1.63633337e-08]\n",
      " [2.32485443e-08 9.99999720e-01 7.10706172e-08 1.85762082e-07]\n",
      " [1.70167600e-08 7.62783780e-08 6.30816344e-01 3.69183563e-01]\n",
      " [2.20645858e-08 2.09297550e-07 3.87559216e-01 6.12440553e-01]]\n",
      "[[9.99999902e-01 3.88515923e-08 2.79467307e-08 3.12758302e-08]\n",
      " [4.64970863e-08 9.99999440e-01 1.87896911e-07 3.25768366e-07]\n",
      " [3.58970937e-08 2.01665219e-07 5.41009751e-01 4.58990011e-01]\n",
      " [4.21728388e-08 3.67042187e-07 4.81835668e-01 5.18163923e-01]]\n",
      "['GAGCACCTCGCCGGTCTGGATCTCTGGCACAGTCAACTAGAGCTATCATG', 'GAGCACCTCTCCGGTCTGGATCGCTGTCACAGTCAACGATAGCTATCAGG']\n",
      "['ACCGATATACGCCATGTAAACGTAGTCTCGGCCACTCTTACCTAAATTACCCTTCGGACAAGCTTCGGAGGACATATACCCCCTACGTCACTGGTGCAAA', 'ACCTATATACTCCATTGAAACTTAGTCTCTTCCACTCGTACCGAAATTACCCTTCGTACAATCTGCTTATTACAGATACCCCCGACGTCACTGTGGCAAA']\n",
      "['GGTACACCACCTAGCGCGCGTTATAAGTCCAAACCAGCGTGACCGACTGTACTACCGAATAACAGCCCGGACGAAGAGAGTATGATTTACAAATACGGAAAAAAAAGCCGCAGTGGCCTAAAGCGGTTGCTGACTCTAAGCTAAATAGTAGGGAGGCGCGAGCTCTTAAAGGATGAATATTATACCAAAACGGTCACACCTCGGGGTTTAAGCCTTATGGGTAGCTTCACCACACGTCCAGACTATTAATCAACTACCAACAACCAGGACATCTGCGTGTGCGCATCGTTCCGACGTGCTCCTGGCGTCCACTTTGTGGTCCTGTTAATGGCCTTTTTTGCTAGACGGACTCCTCATGGTGCCCTAAGATGGAATGGGCAGGCGACCCCTGGAAAGAGAATCAAGTAAATGTCTAGAACGTAATGTGGAGCAGCCGCGACCTTTACATAGCCCATGGAAAGCCACAATAATCACATTTGTTACCAGTAAACACGCGTACCGAGGCCCTAAAAGAGACAACGATGCATTTACGTCTACTCTCCAAATGATCCACAGTGGAGCTAACTGAGTCCACCCAAGAGGATACTGTCCCTTCGCAAATTTACTAAAGACGCGGTGCGAAATAACGGAAACATGTGGAACTTCAGATGCCGTCTCATTGGACATGGGAAGTATGTGGAAAGAGTCAACTATCGAAGGAACAGTACATTTCCAACTTCAAGGAAACAGAAAAACACCAACATCTGACCAATTTGTTCTACTATTGAACTGGGCCCTTGGATGAAATAGGTATTAACACGAGGACTGATGTAGCAAAGGTTCTTAGACTATGCTCACAAAGTGCTAAGGCATCTAAGAACACAGATAGATCGCCTCATCGTGCTGTTTATTTTAGAGAACCCAATATCCCAGCCTCGGTGAGGGCGATCAAGTCGGCTTTGAACCAAACGGCGTCACAAGTACGTTCGGTCAACTCAGCAAAAAGTTATACCGGGCGAAA', 'GGGACACCACCTATCTCTCTTGATAAGGCCAAACCAGCGTGACCTACTGTACTACCTAAGAACAGCCCGTACTAAGATAGTATGAGGGACAAAGACGTAAAAAAAATCCTCAGTGTCCGAAATCGGGTGCTTACTCTAAGCGAAATAGGAGGGATGCTCTAGCGCTTAAATGATGAAGATTAGACCAAAACTGTCACACCGCTGGGGGTAATCCTTAGGGGTAGCGTCACCACACGTCCAGACTAGTAAGCAACTACCAACAACCAGGACAGCGGCGTGTTCGCATCGTGCCGACTGGCGCCTTGCGGCCACGGTGTGGGCCTGTGAATGGCCGTTTTTTCGATACGTACGCCTCATGGTGCCCTAAGATGGAAGGGTCAGGCGACCCCTTGAAATAGAATCAATGAAATTTCTATAACTTAAGGGGTAGCAGCCTCTACCTGTACAGATCCCATGGAAAGCCACAAGAATCACATGTGGTACCATGAAACACGCTTACCTATGCCCGAAAAGATACAACTATTCATGTACGTCTACTCTCCAAAGGATCCACAGGGGAGCTAACTGATTCCACCCAAGATGAGACTTTCCCTTCGCAAAGTGACTAAAGACTCTGGGCTAAAGAACTGAAACAGGGTGAACGTCATATTCCTGCTCATGGGACATGGGAAGTATTGGGAAATAGTCAACGAGCTAAGGAACAGTACATTTCCAACGGCAATGAAACAGAAAAACACCAACAGCGGACCAATGTGGTCTACGATTGAACGGGTCCCGTGGATGAAATATTTATTAACACTAGGACGGATGGATCAAAGGGTCGGAGACGATGCTCACAAAGGTCGAAGTCATCTAAGAACACAGAGAGATCGCCGCATCTGGCGGTTGAGGGTAGATAACCCAATAGCCCAGCCTCGGGGAGTTCTAGCAAGTCGGCTGGTAACCAAACGTCGTCACAAGGACGTTCTGGCAACTCATCAAAAATTTAGACCGGGCGAAA']\n",
      "['CCGAACGACATAAGGCCGACACCACGGCCATGCAATGACGCACTCCACCG', 'CCTAACGACAGAATTCCGACACCACGTCCAGGCAAGTACGCACTCCACCG']\n",
      "['AAGATACGTTAAACGGCAACCATCCCCCCATAAGAATAGTTGGGGCAAAGACTGTGCTGACTCATTACCCATTCCAAGACCTAGCTCATTAGAAGTTTAC', 'AAGAGACTGGAAACGTCAACCATCCCCCCATAAGAATAGGTTGGGCAAAGACTGTGCTGACTCAGGACCCAGGCCAATACCTAGCGCATTATAATGTTAC']\n",
      "['ACAACCGTGGTGGAGGCCGGAGGCACAGCGATAGCAGATTATCTTGGATTAAGGTATGGCAAAAATGTAAGCCATGCCATTTTCTATGTTAACTTAAGGCGTAAATGATCCACTGGTGAGGTGGGTGGAATAATATTGCCTAGGACCAAACAAAGATTGTAGGTCGCTATAGAGGCACTGCGACGTGTACTCCGCTAAACATGGGAGTTTGCTCTAAAATAAACCATTTACAGACTAAGTCAGCCGTCTGCTGGAAAATTAAGATACTCACACCGAACTGGGGAGTATACCTTAAGCATAAATAAATGTGACAATCATGGTCCACCTCAACGTCAGACCTTATTCGTGAAACAGCTATCAATTTGATCCTATCACCAGTAATCATGAAATCAATCGAGTTCGTCGTAGCCGCGAGCTAACCGTACGTAAGGGACAATGTAGATCTACAGCTTAGCAATGGATACCATTTTTAGAACAAACCTTAGACTGGATATGCACGGAAAGATAGGACCCAATACGCAAGTATGACCGGGTCTTGAATCACCGTACACTTTCTCATCTATTGTGCAAGACACACGGGACGAGTGGCTACAGCCACCCTCAGCTATGAAACACAGAACGTCACGGAACGAAGCATGATAACCACACAACTACACGCGATATAGCGAATTATTATAAACGCAGCCAAACACAGTAAAGCGTCAACTGGAATCGACAACATCGTAGAGTCACTGTGCAGGACATCTCTAGTGATTGCTACGATACAACAATCAGTGGGCGAGGTAGTTCCGCTATCTTTTTCTGATGATGGACGCAGACCGCCTGGTGGCTACCCATGTGAACAGATCAGGACCTCTGCAAGACCTCGGTATTTATCTACACGCGCACGACAAGGGATATAGTGGACCATTGATAAATAACAGCGTCGCAATGATCCTGCCTTTAAGTCCATAAACTACCCAAACAGCTTTAAGCGCGTTCCTAAGATGGAACTAGGTAC', 'ACAACCTTGTTGGAGTCCTGAGGCACAGCGATAGCAGATGAGCGTGTAGGAATGTATGGCAAAAATGTAAGCCATGCCAGTTTCTATTTTAACGTAATGCGTAAATGAGCCACGGGTTAGGTTTTGGTAATAAGAGTTCCTAGGACCAAACAAAGATGGTAGTGCGCTATATATGCACGTCGACTGGGACGCCTCTAAACATGGGAGGTTGCTCGAAAATAAACCATGGACATACGAAGGCATCCGTCGGCTGTAAAAGGAATAGACGCACACCGAACGGTGTAGGAGACCTTAATCATAAATAAAGGGGACAATCAGTTTCCACCGCAACGTCAGACCGGAGGCGTGAAACATCTAGCAAGGTTATCCGATCACCAGTAAGCAGTAAATCAATCTATTTCTGCTTAGCCTCTATCTAACCTGACTTAAGTTACAAGGTAGAGCTACAGCGTAGCAATTTAGACCAGTTTTAGAACAAACCTTATACGGTATAGTCACGGAAAGAGATTACCCAATACGCAAGTATGACCTTGTCTGTAATCACCTGACACTGTCTCATCGATGTTGCAAGACACACTGTACGAGGGGCTACAGCCACCCTCATCTATTAAACACATAACTTCACGGAACGAAGCAGGATAACCACACAACGACACGCTATATAGCGAAGTAGTATAAACTCATCCAAACACATGAAAGCGTCAACTTGAATCGACAACAGCTGAGATGCACGGGGCAGGACAGCTCGAGGGATGGCGACTAGACAACAAGCATTTTGCGATGGAGTGCCGCGAGCGGGTTCGTAGGATGTACGCAGACCGCCTGTGGTCTACCCAGGTGAACAGAGCAGTACCTCTGCAATACCTCGGGAGTTAGCTACACGCGCACTACAAGTGATATAGGTGACCATGGAGAAAGAACATCTTCGCAATGATCCGGCCTTTAATTCCAGAAACTACCCAAACATCTGTAAGCTCTTGCCGAATAGGGAACTAGGGAC']\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import random\n",
    "from random import choice\n",
    "import numpy as np\n",
    "import scipy\n",
    "from scipy.optimize import newton\n",
    "from scipy.optimize import brenth\n",
    "from scipy.optimize import root\n",
    "from scipy.optimize import brentq\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.linalg import expm\n",
    "import matplotlib.pyplot as plt \n",
    "import math\n",
    "import rpy2\n",
    "from rpy2.robjects import r\n",
    "import rpy2.robjects.numpy2ri\n",
    "rpy2.robjects.numpy2ri.activate()\n",
    "\n",
    "# Pi = 0-A 1-C  2-G  3-T\n",
    "p = [.298,.249,.232,.221]\n",
    "\n",
    "#Rate = 0AC   1AG    2AT    3CG   4CT   5GT(constrainte)\n",
    "g = [3.821,21.193,2.614,2.567,34.546, 1.000]\n",
    "\n",
    "#human rates\n",
    "#Rate = 0AC   1AG    2AT    3CG   4CT   5GT(constrainte)\n",
    "#g = [2.5*10**-8, 1.6*10**-8, 2.614*10**-8, 2.567*10**-8, 34.546*10**-8, 1.000]\n",
    "\n",
    "nucs = [\"A\",\"C\",\"G\",\"T\"]\n",
    "\"\"\"\n",
    "Q=[A*AA, C*AC, G*AG, T*AT \n",
    "   A*AC, C*CC, G*CG, T*CT\n",
    "   A*AG, C*GC, G*GG, T*GT\n",
    "   A*AT, C*TC, G*TG, T*TT]\n",
    "\"\"\"\n",
    "AA = float(-(p[1]*g[0]+p[2]*g[1]+p[3]*g[2]))\n",
    "CC = float(-(p[0]*g[0]+p[2]*g[3]+p[3]*g[4]))\n",
    "GG = float(-(p[0]*g[1]+p[1]*g[3]+p[3]*g[5]))\n",
    "TT = float(-(p[0]*g[2]+p[1]*g[4]+p[2]*g[5]))\n",
    "\n",
    "A = np.array(\n",
    "     [[AA,p[1]*g[0],p[2]*g[1],p[3]*g[2]],\n",
    "     [p[0]*g[0],CC,p[2]*g[3],p[3]*g[4]],\n",
    "     [p[0]*g[1],p[1]*g[3],GG,p[3]*g[5]],\n",
    "     [p[0]*g[2],p[1]*g[4],p[2]*g[5], TT]]\n",
    ")\n",
    "\n",
    "# A_s = -divide(A,(.25*(-1/16.0 * 3))*4)\n",
    "# A_s\n",
    "\n",
    "# print(A)\n",
    "\n",
    "norm_factor = -(p[0]*AA + p[1]*CC + p[2]*GG + p[3]*TT)\n",
    "\n",
    "A_s = (np.divide(A, norm_factor))\n",
    "\n",
    "print(\"As=\", A_s)\n",
    "\n",
    "P = expm(np.multiply(A_s, 0.01))\n",
    "\n",
    "def RandomSeq(length, p):\n",
    "    DNA = \"\"\n",
    "    for count in range(length):\n",
    "        DNA+=np.random.choice(nucs, p = p)\n",
    "        #DNA+=np.random.choice(nucs, p) \n",
    "    return DNA\n",
    "\n",
    "rand50 = RandomSeq(50,p)\n",
    "# print(rand50)\n",
    "\n",
    "\n",
    "def RandomAlignment(length, p, P):\n",
    "    \n",
    "    DNA1 = \"\"\n",
    "    DNA2 = \"\"\n",
    "    for count in range(length):\n",
    "        nuc1 = np.random.choice(nucs, p = p)\n",
    "        DNA1+=nuc1\n",
    "        #DNA+=np.random.choice(nucs, p)\n",
    "        if nuc1 == \"A\":\n",
    "            nuc2 = np.random.choice(nucs, p = P[0])\n",
    "        elif nuc1 == \"C\":\n",
    "            nuc2 = np.random.choice(nucs, p = P[1])\n",
    "        elif nuc1 == \"G\":\n",
    "            nuc2 = np.random.choice(nucs, p = P[2])\n",
    "        elif nuc1 == \"T\":\n",
    "            nuc2 = np.random.choice(nucs, p = P[3])\n",
    "        else:\n",
    "            RaiseError(\"ERROR\")\n",
    "       \n",
    "        DNA2+=nuc2\n",
    "\n",
    "    alignment = [DNA1, DNA2]\n",
    "    return alignment;\n",
    "\n",
    "P1 = np.linalg.matrix_power(P, 1)\n",
    "print(\"P =\", P)\n",
    "P10 = np.linalg.matrix_power(P, 10)\n",
    "P1024 = np.linalg.matrix_power(P, 1024)\n",
    "\n",
    "# Usage: RandomAlignment(Length of alignmemt, stationary distribution, P matrix, n)\n",
    "# RandomAlignment(50, p, np.linalg.matrix_power(P, 16))\n",
    "\n",
    "#set your n\n",
    "n = 32\n",
    "print(np.linalg.matrix_power(P, n))\n",
    "print(np.linalg.matrix_power(P, 2*n))\n",
    "\n",
    "P_matrix = np.linalg.matrix_power(P, n)\n",
    "\n",
    "randTest50 = RandomAlignment(50, p, np.linalg.matrix_power(P, n))\n",
    "print(randTest50)\n",
    "\n",
    "randTest100 = RandomAlignment(100, p, np.linalg.matrix_power(P, n))\n",
    "print(randTest100)\n",
    "\n",
    "randTest1000 = RandomAlignment(1000, p, np.linalg.matrix_power(P, n))\n",
    "print(randTest1000)\n",
    "\n",
    "n2randTest50 = RandomAlignment(50, p, np.linalg.matrix_power(P, 2*n))\n",
    "print(n2randTest50)\n",
    "\n",
    "n2randTest100 = RandomAlignment(100, p, np.linalg.matrix_power(P, 2*n))\n",
    "print(n2randTest100)\n",
    "\n",
    "n2randTest1000 = RandomAlignment(1000, p, np.linalg.matrix_power(P, 2*n))\n",
    "print(n2randTest1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  121. -1637. -1668. -1642.]\n",
      " [-1637.   139. -1500. -1399.]\n",
      " [-1668. -1500.   100.    51.]\n",
      " [-1642. -1399.    51.   102.]]\n",
      "[[  1.21066174 -16.36636152 -16.6784053  -16.41863017]\n",
      " [-16.36636152   1.3903021  -14.99857394 -13.98920653]\n",
      " [-16.6784053  -14.99857394   1.00027739   0.51313128]\n",
      " [-16.41863017 -13.98920653   0.51313128   1.01928918]]\n",
      "[1.e-05 2.e+00]\n",
      "Diag(p) [[0.298 0.    0.    0.   ]\n",
      " [0.    0.249 0.    0.   ]\n",
      " [0.    0.    0.232 0.   ]\n",
      " [0.    0.    0.    0.221]]\n",
      "SD exp*3.005106 [[8.26447658e+157 0.00000000e+000 0.00000000e+000 0.00000000e+000]\n",
      " [0.00000000e+000 2.56467603e+181 0.00000000e+000 0.00000000e+000]\n",
      " [0.00000000e+000 0.00000000e+000 3.23664698e+130 3.63202325e+066]\n",
      " [0.00000000e+000 0.00000000e+000 3.63202325e+066 1.31915929e+133]]\n",
      "[[0.088804 0.       0.       0.      ]\n",
      " [0.       0.062001 0.       0.      ]\n",
      " [0.       0.       0.053824 0.      ]\n",
      " [0.       0.       0.       0.048841]]\n",
      "[[-6.61573540e-04 -9.99999977e-01 -9.99999983e-01 -9.99999978e-01]\n",
      " [-9.99999981e-01 -3.02336804e-04 -9.99999924e-01 -9.99999791e-01]\n",
      " [-9.99999987e-01 -9.99999929e-01 -3.69358616e-01 -6.13652443e-01]\n",
      " [-9.99999984e-01 -9.99999814e-01 -6.31970646e-01 -3.87123957e-01]]\n",
      "[[  363.617826 -4919.358522 -5012.516808 -4934.384052]\n",
      " [-4919.358522   417.709734 -4507.659    -4204.143294]\n",
      " [-5012.516808 -4507.659      300.5106     153.260406]\n",
      " [-4934.384052 -4204.143294   153.260406   306.520812]]\n",
      "[[3.35348465e+00 7.77319086e-08 5.70122318e-08 7.39408787e-08]\n",
      " [7.77319086e-08 4.01485005e+00 3.05902321e-07 8.39885722e-07]\n",
      " [5.70122318e-08 3.05902321e-07 2.71828183e+00 1.66529119e+00]\n",
      " [7.39408787e-08 8.39885722e-07 1.66529119e+00 2.77319476e+00]]\n",
      "lambda= [1] 0.009986909\n",
      "\n"
     ]
    }
   ],
   "source": [
    "S = np.log(np.matmul(P_matrix, np.diag(np.divide(1,p))))\n",
    "Sd = np.round(100*S)\n",
    "print(Sd)\n",
    "print(S)\n",
    "t = np.array([0.00001, 2])\n",
    "print(t)\n",
    "p = p\n",
    "print(\"Diag(p)\", np.diag(p))\n",
    "print(\"SD exp*3.005106\", np.exp(3.005106*Sd))\n",
    "\n",
    "print(np.matmul(np.diag(p), np.diag(p), np.diag(p)))\n",
    "print((np.matmul(np.diag(p), np.exp(0.01*Sd), np.diag(p)))-1)\n",
    "\n",
    "print(3.005106*Sd)\n",
    "print(np.exp(0.01*Sd))\n",
    "#print(np.sum(np.matmul(np.diag(p), ))\n",
    "\n",
    "#for i in np.arange(0.0001, 2, 0.1):\n",
    "#    print(i)\n",
    "#for i in range(0, 10):\n",
    "#    print(np.sum((np.matmul(np.diag(p), np.exp(i*S), np.diag(p))))-1)\n",
    "\n",
    "#def mgf(t, p, S):\n",
    "#    m = np.sum((np.matmul(np.diag(p), np.exp(t*S), np.diag(p))))-1\n",
    "#    return m;\n",
    "\n",
    "rmgf = r(\"function(t,p,S) {sum(diag(p) %*% exp(t*S) %*% diag(p)) - 1}\")\n",
    "\n",
    "lamb = r.uniroot(rmgf, interval=t, p=p, S=Sd)[0]\n",
    "\n",
    "print(\"lambda=\",lamb)\n",
    "\n",
    "#v_mgf = np.vectorize(mgf, otypes=[float])\n",
    "#print(v_mgf)\n",
    "#\n",
    "#args=(p, Sd)\n",
    "#print(args)\n",
    "#\n",
    "#l = fsolve(v_mgf, [0.001, 2], args=args)\n",
    "#l = fsolve(lambda x: v_mgf(x, p, Sd), [0.001, 2])\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.009986909345850475\n",
      "0.64354\n",
      "4.510461294965931e-06\n",
      "0.6435370973377382 0.009937205766037072 107.8122865646166\n",
      "1.0713528756995492\n",
      "C= 0.6006770616241032\n",
      "A= 1.5359428571401523\n"
     ]
    }
   ],
   "source": [
    "S_vector = Sd.flatten()\n",
    "#print(S_vector)\n",
    "pi = np.array([[p]])\n",
    "pi2 = np.transpose(pi)\n",
    "#print(pi)\n",
    "#print(pi2)\n",
    "\n",
    "null_probs = (pi*pi2)\n",
    "null_probs = null_probs.flatten()\n",
    "#print(\"null probs=\", null_probs)\n",
    "#print(\"Score vector=\", S_vector)\n",
    "reps = 100000\n",
    "H = np.repeat(0, reps)\n",
    "L = np.repeat(-1, reps)\n",
    "Qt = np.repeat(0, max(S_vector)+1)\n",
    "R = np.repeat(0, -(min(S_vector)-1))\n",
    "#print(R, L, H)\n",
    "exc_score = 0\n",
    "i = 1\n",
    "while i < reps:\n",
    "    exc_score = exc_score + np.random.choice(S_vector, p=null_probs)\n",
    "    L[i] = L[i]+1\n",
    "    if (exc_score > 0) and (L[i] == 0):\n",
    "        Qt[int(exc_score)] = Qt[int(exc_score)]+1\n",
    "    if H[i] < exc_score:\n",
    "        H[i] = exc_score\n",
    "    if exc_score < 0:\n",
    "        R[int(-exc_score)] = R[int(-exc_score)]+1\n",
    "        i = i + 1\n",
    "        exc_score = 0\n",
    "\n",
    "        \n",
    "R_total = list(R/reps)\n",
    "#print(R_total)\n",
    "print(lamb[0])\n",
    "Qt = list(Qt/reps)\n",
    "Q_avg = float(1-np.sum(Qt))\n",
    "#print(type(Q_avg))\n",
    "print(Q_avg)\n",
    "#calculate C as shown in notes\n",
    "C1_a = 0\n",
    "for i in range(1, len(R_total)):\n",
    "    C1_a += (R_total[i-1]*(np.exp(-(lamb[0])*i)))\n",
    "    #print(C1_a)\n",
    "\n",
    "print(C1_a)\n",
    "C1 = Q_avg*(1-C1_a)    \n",
    "C2 = (1-np.exp(-lamb[0]))\n",
    "C3 = 0\n",
    "for i in range(1, len(Qt)):\n",
    "    C3 += i*Qt[i]*np.exp(i*lamb[0])\n",
    "    \n",
    "    \n",
    "print(C1, C2, C3)\n",
    "\n",
    "C = (C1/(C2*C3))\n",
    "print(C2*C3)\n",
    "\n",
    "print(\"C=\", C)\n",
    "\n",
    "A1 = 0\n",
    "for i in range(1, len(R)):\n",
    "    A1 += i*R_total[i]\n",
    "A1 = -(A1)\n",
    "\n",
    "A2 = 0\n",
    "for i in range(1, len(S_vector)):\n",
    "    A2 += S_vector[i]*null_probs[i]\n",
    "\n",
    "A = A1/A2\n",
    "print(\"A=\", A)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, -8, -7, -17, -16, -26, -25, -35, -34, -44, -43, -42, -41, -40, -39, -49, -48, -47, -46, -45, -44, -43, -42, -41, -40, -39, -38, -37, -36, -35, -34, -33, -43, -42, -41, -40, -39, -38, -37, -36, -35, -34, -33, -32, -42, -41, -40, -50, -49, -48, -47, -46, -45, -44, -43, -42, -41, -51, -50, -49, -59, -58, -57, -56, -55, -65, -64, -63, -62, -61, -60, -59, -58, -68, -78, -88, -87, -86, -85, -84, -83, -93, -92, -91, -90, -100, -99, -98, -97, -96, -95, -94, -93, -92, -102, -101, -100, -110, -109, -108, -107, -106, -105, -115, -114, -113, -123, -122, -121, -120, -130, -129, -128, -127, -137, -136, -135, -134, -133, -143, -142, -141, -140, -139, -138, -137, -136, -135, -134, -144, -143, -142, -141, -140, -139, -138, -148, -147, -146, -145, -144, -143, -153, -152, -151, -161, -160, -170, -169, -168, -167, -177, -176, -175, -174, -173, -172, -171, -181, -180, -179, -178, -177, -176, -175, -185, -184, -183, -182, -181, -191, -190, -189, -188, -187, -186, -185, -184, -183, -193, -192, -191, -190, -189, -188, -187, -186, -185, -195, -194, -204, -203, -202, -201, -211, -221, -220, -219, -218, -228, -227, -226, -225, -224, -223, -233, -232, -231, -230, -229, -228, -227, -226, -236, -235, -234, -233, -232, -231, -230, -229, -228, -227, -226, -225, -224, -223, -222, -221, -220, -219, -218, -217, -227, -226, -225, -224, -234, -233, -232, -231, -230, -229, -228, -227, -226, -225, -224, -223, -222, -221, -220, -219, -218, -217, -216, -215, -214, -213, -223, -222, -232, -231, -230, -229, -228, -227, -226, -236, -235, -234, -233, -232, -231, -230, -229, -228, -238, -237, -236, -235, -234, -233, -243, -253, -252, -251, -261, -260, -259, -258, -268, -267, -266, -265, -275, -274, -273, -272, -271, -281, -291, -290, -289, -288, -287, -286, -296, -295, -294, -293, -292, -291, -301, -300, -299, -298, -297, -296, -295, -294, -304, -303, -302, -301, -300, -299, -309, -308, -318, -317, -327, -326, -325, -324, -334, -333, -332, -342, -341, -340, -339, -338, -337, -336, -335, -334, -333, -332, -331, -330, -329, -328, -327, -326, -325, -324, -323, -322, -321, -320, -319, -329, -328, -327, -337, -336, -335, -334, -333, -332, -331, -330, -329, -328, -327, -326, -325, -335, -334, -333, -332, -331, -341, -340, -339, -338, -337, -336, -335, -334, -333, -343, -353, -352, -351, -350, -349, -359, -358, -357, -356, -355, -365, -364, -363, -362, -372, -371, -370, -369, -379, -378, -388, -387, -397, -396, -395, -394, -393, -392, -391, -390, -400, -399, -409, -408, -407, -406, -405, -415, -414, -413, -412, -411, -421, -420, -430, -429, -428, -427, -426, -425, -424, -423, -422, -421, -420, -419, -418, -417, -416, -415, -414, -413, -423, -422, -421, -420, -419, -418, -417, -416, -415, -425, -424, -423, -433, -432, -431, -430, -429, -428, -438, -448, -447, -446, -445, -444, -443, -442, -441, -440, -450, -449, -448, -447, -446, -456, -455, -465, -464, -463, -462, -461, -471, -470, -469, -468, -467, -466, -465, -475, -474, -473, -472, -471, -470, -480, -479, -478, -488, -487, -486, -485, -495, -494, -493, -492, -491, -490, -489, -488, -487, -486, -485, -484, -483, -482, -481, -480, -479, -478, -488, -487, -486, -485, -484, -483, -482, -481, -480, -479, -489, -488, -487, -486, -485, -484, -483, -482, -481, -480, -479, -478, -477, -487, -486, -485, -484, -483, -482, -481, -480, -479, -478, -477, -476, -486, -485, -484, -494, -493, -492, -491, -501, -500, -499, -498, -497, -496, -495, -494, -493, -492, -491, -490, -489, -499, -498, -508, -507, -506, -505, -504, -503, -502, -501, -500, -499, -509, -508, -518, -517, -527, -526, -525, -535, -534, -533, -532, -542, -541, -540, -539, -549, -548, -547, -546, -545, -544, -543, -553, -552, -562, -572, -571, -570, -569, -568, -578, -577, -576, -575, -585, -584, -583, -593, -592, -591, -601, -611, -610, -609, -608, -607, -606, -616, -615, -614, -613, -612, -611, -610, -609, -608, -607, -606, -605, -604, -603, -602, -601, -611, -621, -620, -619, -618, -617, -616, -626, -625, -624, -623, -622, -621, -620, -619, -629, -628, -638, -637, -647, -646, -645, -644, -643, -642, -641, -640, -639, -638, -637, -636, -635, -634, -633, -632, -631, -630, -629, -628, -627, -626, -636, -646, -645, -644, -643, -653, -652, -651, -650, -649, -648, -647, -646, -645, -644, -643, -642, -641, -640, -639, -638, -637, -636, -635, -634, -633, -643, -642, -652, -651, -650, -649, -648, -647, -646, -645, -655, -654, -653, -663, -662, -661, -660, -659, -658, -668, -667, -666, -665, -664, -663, -662, -661, -671, -670, -669, -679, -678, -677, -676, -686, -685, -684, -683, -682, -681, -680, -679, -678, -677, -676, -675, -685, -695, -694, -693, -692, -691, -690, -689, -688, -687, -686, -696, -695, -694, -693, -692, -691, -701, -700, -699, -698, -697, -707, -706, -716, -715, -714, -713, -712, -711, -710, -720, -719, -718, -728, -738, -737, -736, -735, -734, -744, -743, -742, -741, -740, -739, -738, -737, -736, -735, -734, -733, -732, -742, -752, -751, -761, -760, -759, -758, -768, -767, -766, -765, -764, -763, -762, -761, -760, -759, -758, -757, -756, -755, -754, -753, -752, -762, -761, -760, -759, -758, -757, -756, -755, -754, -764, -763, -762, -761, -760, -770, -780, -779, -778, -788, -787, -786, -785, -795, -794, -804, -814, -824, -823, -822, -821, -820, -830, -829, -828, -827, -826, -825, -824, -823, -822, -821, -831, -830, -829, -828, -827, -826, -825, -824, -823, -822, -821, -820, -830, -829, -828, -827, -837, -847, -846, -856, -855, -865, -864, -863, -862, -861, -860, -859, -858, -857, -856, -855, -865, -875, -885, -884, -883, -882, -881, -880, -879, -878, -877, -876, -886, -885, -884, -883, -882, -881, -880, -879, -878, -877, -887, -886, -885, -884, -883, -882, -881, -891, -890, -900, -899, -898, -897, -896, -895, -894, -893, -903, -902, -901, -900, -899, -898, -897, -907, -906, -905, -904, -914, -913, -912, -911, -910, -909, -908, -907, -906, -905, -904, -903]\n",
      "[1, 2, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 0, 1, 2, 0, 1, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 0, 1, 2, 3, 0, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 9, 10, 11, 12, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 13, 14, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 0, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 4, 0, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 0, 1, 0, 1, 0, 1, 2, 3, 0, 1, 2, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 13, 14, 15, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 7, 8, 9, 10, 11, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 0, 1, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 7, 8, 9, 10, 11, 12, 13, 14, 15, 5, 6, 7, 0, 1, 2, 3, 4, 5, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 0, 1, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 9, 10, 11, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 2, 3, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 0, 1, 0, 0, 1, 2, 3, 4, 0, 1, 2, 3, 0, 1, 2, 0, 1, 2, 0, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 0, 1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 11, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 10, 11, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 0, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 2, 0, 1, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 6, 7, 8, 9, 10, 11, 12, 13, 14, 4, 5, 6, 7, 8, 0, 0, 1, 2, 0, 1, 2, 3, 0, 1, 0, 0, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 1, 2, 3, 4, 0, 0, 1, 0, 1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 0, 1, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]\n",
      "['N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'Y', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'Y', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'Y', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'Y', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'Y', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'Y', 'Y', 'Y', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'Y', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N']\n",
      "Match, no match = 827 173\n",
      "23\n"
     ]
    }
   ],
   "source": [
    "def AlignmentScorer(list, match_score, no_match_score):\n",
    "    DNA1 = list[0]\n",
    "    DNA2 = list[1]\n",
    "    match = 0\n",
    "    no_match = 0\n",
    "    score = []\n",
    "    excursion = []\n",
    "    ladder = []\n",
    "    max_height = []\n",
    "    low = 0\n",
    "    for i in range(0, len(DNA1)):\n",
    "        if DNA1[i] == DNA2[i]:\n",
    "            #print(DNA1[i], DNA2[i], \"MATCH\")\n",
    "            #print(i)\n",
    "            match += 1  \n",
    "            ladder.append(\"N\")\n",
    "            if i > 0:\n",
    "                score.append(match_score+score[i-1])\n",
    "            elif i == 0:\n",
    "                score.append(match_score + low)\n",
    "            \n",
    "            excursion.append(abs(low-score[i]))\n",
    "            \n",
    "        else:\n",
    "            #print(DNA1[i], DNA2[i], \"NO MATCH\")\n",
    "            no_match+=1\n",
    "            if i > 0:\n",
    "                score.append(no_match_score+score[i-1])\n",
    "            elif i == 0:\n",
    "                score.append(no_match_score+low)\n",
    "            if score[i] < low:\n",
    "                #print(score[i])\n",
    "                ladder.append(\"Y\")\n",
    "                low = score[i]\n",
    "                excursion.append(abs(low-score[i]))\n",
    "            else:\n",
    "                ladder.append(\"N\")\n",
    "                excursion.append(abs(low-score[i]))\n",
    "    \n",
    "    alignment_results = [score, excursion, ladder, match, no_match]\n",
    "    \n",
    "    return alignment_results;\n",
    "    \n",
    "y = AlignmentScorer(randTest1000, 1, -10)\n",
    "print(y[0])\n",
    "print(y[1])\n",
    "print(y[2])\n",
    "print(\"Match, no match =\", y[3], y[4])\n",
    "#print(P)\n",
    "\n",
    "max_height= np.max(y[1])\n",
    "length = len(y[1])\n",
    "print(max_height)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of sequence = 10000\n",
      "Max Height =  89\n",
      "E= 1607.8661152524585\n"
     ]
    }
   ],
   "source": [
    "l = 10000\n",
    "n = 10000\n",
    "P_matrix = np.linalg.matrix_power(P, n)\n",
    "x = RandomAlignment(l, p, P_matrix)\n",
    "#print(x[0])\n",
    "#print(x[1])\n",
    "y = AlignmentScorer(x, 1, -4)\n",
    "max_height= np.max(y[1])\n",
    "length = len(y[1])\n",
    "E = length*(C*np.exp(-lamb[0]*max_height))/A\n",
    "print(\"Length of sequence =\", length)\n",
    "print(\"Max Height = \", max_height)\n",
    "print(\"E=\", E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}