tests/phylogenetic_invariants.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Phylogenetic Invariants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AAAA AAAC AAAG AAAT  AACA AACC AACG AACT  AAGA AAGC AAGG AAGT  AATA AATC AATG AATT\n",
      "ACAA ACAC ACAG ACAT  ACCA ACCC ACCG ACCT  ACGA ACGC ACGG ACGT  ACTA ACTC ACTG ACTT\n",
      "AGAA AGAC AGAG AGAT  AGCA AGCC AGCG AGCT  AGGA AGGC AGGG AGGT  AGTA AGTC AGTG AGTT\n",
      "ATAA ATAC ATAG ATAT  ATCA ATCC ATCG ATCT  ATGA ATGC ATGG ATGT  ATTA ATTC ATTG ATTT\n",
      "\n",
      "CAAA CAAC CAAG CAAT  CACA CACC CACG CACT  CAGA CAGC CAGG CAGT  CATA CATC CATG CATT\n",
      "CCAA CCAC CCAG CCAT  CCCA CCCC CCCG CCCT  CCGA CCGC CCGG CCGT  CCTA CCTC CCTG CCTT\n",
      "CGAA CGAC CGAG CGAT  CGCA CGCC CGCG CGCT  CGGA CGGC CGGG CGGT  CGTA CGTC CGTG CGTT\n",
      "CTAA CTAC CTAG CTAT  CTCA CTCC CTCG CTCT  CTGA CTGC CTGG CTGT  CTTA CTTC CTTG CTTT\n",
      "\n",
      "GAAA GAAC GAAG GAAT  GACA GACC GACG GACT  GAGA GAGC GAGG GAGT  GATA GATC GATG GATT\n",
      "GCAA GCAC GCAG GCAT  GCCA GCCC GCCG GCCT  GCGA GCGC GCGG GCGT  GCTA GCTC GCTG GCTT\n",
      "GGAA GGAC GGAG GGAT  GGCA GGCC GGCG GGCT  GGGA GGGC GGGG GGGT  GGTA GGTC GGTG GGTT\n",
      "GTAA GTAC GTAG GTAT  GTCA GTCC GTCG GTCT  GTGA GTGC GTGG GTGT  GTTA GTTC GTTG GTTT\n",
      "\n",
      "TAAA TAAC TAAG TAAT  TACA TACC TACG TACT  TAGA TAGC TAGG TAGT  TATA TATC TATG TATT\n",
      "TCAA TCAC TCAG TCAT  TCCA TCCC TCCG TCCT  TCGA TCGC TCGG TCGT  TCTA TCTC TCTG TCTT\n",
      "TGAA TGAC TGAG TGAT  TGCA TGCC TGCG TGCT  TGGA TGGC TGGG TGGT  TGTA TGTC TGTG TGTT\n",
      "TTAA TTAC TTAG TTAT  TTCA TTCC TTCG TTCT  TTGA TTGC TTGG TTGT  TTTA TTTC TTTG TTTT\n",
      "\n",
      "[['AAAA' 'AAAC' 'AAAG' 'AAAT']\n",
      " ['ACAA' 'ACAC' 'ACAG' 'ACAT']\n",
      " ['AGAA' 'AGAC' 'AGAG' 'AGAT']\n",
      " ['ATAA' 'ATAC' 'ATAG' 'ATAT']]\n"
     ]
    }
   ],
   "source": [
    "pi_string = \"\"\"\\\n",
    "AAAA AAAC AAAG AAAT  AACA AACC AACG AACT  AAGA AAGC AAGG AAGT  AATA AATC AATG AATT\n",
    "ACAA ACAC ACAG ACAT  ACCA ACCC ACCG ACCT  ACGA ACGC ACGG ACGT  ACTA ACTC ACTG ACTT\n",
    "AGAA AGAC AGAG AGAT  AGCA AGCC AGCG AGCT  AGGA AGGC AGGG AGGT  AGTA AGTC AGTG AGTT\n",
    "ATAA ATAC ATAG ATAT  ATCA ATCC ATCG ATCT  ATGA ATGC ATGG ATGT  ATTA ATTC ATTG ATTT\n",
    "\n",
    "CAAA CAAC CAAG CAAT  CACA CACC CACG CACT  CAGA CAGC CAGG CAGT  CATA CATC CATG CATT\n",
    "CCAA CCAC CCAG CCAT  CCCA CCCC CCCG CCCT  CCGA CCGC CCGG CCGT  CCTA CCTC CCTG CCTT\n",
    "CGAA CGAC CGAG CGAT  CGCA CGCC CGCG CGCT  CGGA CGGC CGGG CGGT  CGTA CGTC CGTG CGTT\n",
    "CTAA CTAC CTAG CTAT  CTCA CTCC CTCG CTCT  CTGA CTGC CTGG CTGT  CTTA CTTC CTTG CTTT\n",
    "\n",
    "GAAA GAAC GAAG GAAT  GACA GACC GACG GACT  GAGA GAGC GAGG GAGT  GATA GATC GATG GATT\n",
    "GCAA GCAC GCAG GCAT  GCCA GCCC GCCG GCCT  GCGA GCGC GCGG GCGT  GCTA GCTC GCTG GCTT\n",
    "GGAA GGAC GGAG GGAT  GGCA GGCC GGCG GGCT  GGGA GGGC GGGG GGGT  GGTA GGTC GGTG GGTT\n",
    "GTAA GTAC GTAG GTAT  GTCA GTCC GTCG GTCT  GTGA GTGC GTGG GTGT  GTTA GTTC GTTG GTTT\n",
    "\n",
    "TAAA TAAC TAAG TAAT  TACA TACC TACG TACT  TAGA TAGC TAGG TAGT  TATA TATC TATG TATT\n",
    "TCAA TCAC TCAG TCAT  TCCA TCCC TCCG TCCT  TCGA TCGC TCGG TCGT  TCTA TCTC TCTG TCTT\n",
    "TGAA TGAC TGAG TGAT  TGCA TGCC TGCG TGCT  TGGA TGGC TGGG TGGT  TGTA TGTC TGTG TGTT\n",
    "TTAA TTAC TTAG TTAT  TTCA TTCC TTCG TTCT  TTGA TTGC TTGG TTGT  TTTA TTTC TTTG TTTT\n",
    "\"\"\"\n",
    "\n",
    "## print as string\n",
    "print pi_string\n",
    "\n",
    "## also save as indexable array\n",
    "pi_arr = np.array(pi_string.split()).reshape(16,16)\n",
    "print pi_arr[:4,:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a simple alignment for four taxa\n",
    "Ordered so the first two rows have a more similar sequence than the latter two (e.g., see CATCAT vs. GATGAT). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['C' 'A' 'T' 'C' 'A' 'T' 'C' 'A' 'T' 'C']\n",
      " ['C' 'A' 'T' 'C' 'A' 'T' 'G' 'A' 'T' 'C']\n",
      " ['G' 'A' 'T' 'G' 'A' 'T' 'G' 'A' 'T' 'C']\n",
      " ['G' 'A' 'T' 'G' 'A' 'T' 'G' 'A' 'T' 'C']] ...\n"
     ]
    }
   ],
   "source": [
    "seqalign = np.array([\n",
    "    list(\"CATCATCATCAT\")*10,\n",
    "    list(\"CATCATGATCAT\")*10,\n",
    "    list(\"GATGATGATCCC\")*10,\n",
    "    list(\"GATGATGATCCC\")*10])\n",
    "\n",
    "print seqalign[:, :10], '...'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A trick for dealing with large seq arrays\n",
    "We convert it to integers because they are easier to deal with than characters. The 'view' argument is used here to convert chars to 8bit integers. I then reassign the ints as 0, 1, 2, or 3, which can be used to quickly fill the the count array using a set of rules designated below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['C' 'A' 'T' 'C' 'A' 'T' 'C' 'A' 'T' 'C']\n",
      " ['C' 'A' 'T' 'C' 'A' 'T' 'G' 'A' 'T' 'C']\n",
      " ['G' 'A' 'T' 'G' 'A' 'T' 'G' 'A' 'T' 'C']\n",
      " ['G' 'A' 'T' 'G' 'A' 'T' 'G' 'A' 'T' 'C']]\n",
      "[[67 65 84 67 65 84 67 65 84 67]\n",
      " [67 65 84 67 65 84 71 65 84 67]\n",
      " [71 65 84 71 65 84 71 65 84 67]\n",
      " [71 65 84 71 65 84 71 65 84 67]]\n",
      "[[1 0 3 1 0 3 1 0 3 1]\n",
      " [1 0 3 1 0 3 2 0 3 1]\n",
      " [2 0 3 2 0 3 2 0 3 1]\n",
      " [2 0 3 2 0 3 2 0 3 1]]\n"
     ]
    }
   ],
   "source": [
    "## print as char\n",
    "print seqalign[:, :10]\n",
    "\n",
    "## get a copy with dtype=int8\n",
    "seqint = np.copy(seqalign.view(np.int8))\n",
    "print seqint[:, :10]\n",
    "\n",
    "## convert ints to indexes of matrix\n",
    "seqint[seqint == 65] = 0    ## As\n",
    "seqint[seqint == 67] = 1    ## Cs\n",
    "seqint[seqint == 71] = 2    ## Ts\n",
    "seqint[seqint == 84] = 3    ## Gs\n",
    "\n",
    "print seqint[:, :10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fill the invariants counts array\n",
    "This records the number of times each of the patterns in this array are observed. We will also fill two \"alternate arrays\", representing the site patterns observed for other two possible resolutions of our quartet. This first array assumes taxa are ordered as 12|34. We can fill the two alternates (13|24 and 14|23) from the patterns observed here. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[30  0  0  0  0 10  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0 10  0  0  0  0 20  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0 10  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0 10  0  0  0  0  0  0  0  0  0 30]]\n"
     ]
    }
   ],
   "source": [
    "## create a 16x16 array of zeros\n",
    "invcounts = np.zeros((16,16), dtype=int)\n",
    "\n",
    "## fill the array according to a set of rules\n",
    "for idx in xrange(seqalign.shape[1]):\n",
    "    i = seqint[:, idx] \n",
    "    invcounts[(4*i[0])+i[1], (4*i[2])+i[3]] += 1\n",
    "    \n",
    "print invcounts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fill the three alternative matrices with site counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[30  0  0  0  0 10  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0 10  0  0  0  0 20  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0 10  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0 10  0  0  0  0  0  0  0  0  0 30]]\n",
      "\n",
      " [[30  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0 10  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0 10  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0 20  0  0  0 10  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0 10  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 30]]\n",
      "\n",
      " [[30  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0 10  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0 10  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0 20  0  0  0 10  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0 10  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      "  [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 30]]]\n"
     ]
    }
   ],
   "source": [
    "## Let's create three arrays \n",
    "mats_ints = np.zeros((3, 16, 16), dtype=int)\n",
    "mats_ints[0] = invcounts\n",
    "\n",
    "## Fill the alternates\n",
    "x = 0\n",
    "for i in xrange(0, 16, 4):\n",
    "    for j in xrange(0, 16, 4):\n",
    "        mats_ints[1, i:i+4, j:j+4] = mats_ints[0, x].reshape(4, 4)\n",
    "        mats_ints[2, i:i+4, j:j+4] = mats_ints[0, x].reshape(4, 4).T\n",
    "        x += 1\n",
    "\n",
    "print mats_ints[:, :, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What does this look like as strings?\n",
    "We print this just for our edification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[['AAAA' 'AAAC' 'AAAG' 'AAAT' 'AACA' 'AACC' 'AACG' 'AACT' 'AAGA' 'AAGC']\n",
      "  ['ACAA' 'ACAC' 'ACAG' 'ACAT' 'ACCA' 'ACCC' 'ACCG' 'ACCT' 'ACGA' 'ACGC']\n",
      "  ['AGAA' 'AGAC' 'AGAG' 'AGAT' 'AGCA' 'AGCC' 'AGCG' 'AGCT' 'AGGA' 'AGGC']\n",
      "  ['ATAA' 'ATAC' 'ATAG' 'ATAT' 'ATCA' 'ATCC' 'ATCG' 'ATCT' 'ATGA' 'ATGC']\n",
      "  ['CAAA' 'CAAC' 'CAAG' 'CAAT' 'CACA' 'CACC' 'CACG' 'CACT' 'CAGA' 'CAGC']\n",
      "  ['CCAA' 'CCAC' 'CCAG' 'CCAT' 'CCCA' 'CCCC' 'CCCG' 'CCCT' 'CCGA' 'CCGC']\n",
      "  ['CGAA' 'CGAC' 'CGAG' 'CGAT' 'CGCA' 'CGCC' 'CGCG' 'CGCT' 'CGGA' 'CGGC']\n",
      "  ['CTAA' 'CTAC' 'CTAG' 'CTAT' 'CTCA' 'CTCC' 'CTCG' 'CTCT' 'CTGA' 'CTGC']\n",
      "  ['GAAA' 'GAAC' 'GAAG' 'GAAT' 'GACA' 'GACC' 'GACG' 'GACT' 'GAGA' 'GAGC']\n",
      "  ['GCAA' 'GCAC' 'GCAG' 'GCAT' 'GCCA' 'GCCC' 'GCCG' 'GCCT' 'GCGA' 'GCGC']\n",
      "  ['GGAA' 'GGAC' 'GGAG' 'GGAT' 'GGCA' 'GGCC' 'GGCG' 'GGCT' 'GGGA' 'GGGC']\n",
      "  ['GTAA' 'GTAC' 'GTAG' 'GTAT' 'GTCA' 'GTCC' 'GTCG' 'GTCT' 'GTGA' 'GTGC']\n",
      "  ['TAAA' 'TAAC' 'TAAG' 'TAAT' 'TACA' 'TACC' 'TACG' 'TACT' 'TAGA' 'TAGC']\n",
      "  ['TCAA' 'TCAC' 'TCAG' 'TCAT' 'TCCA' 'TCCC' 'TCCG' 'TCCT' 'TCGA' 'TCGC']\n",
      "  ['TGAA' 'TGAC' 'TGAG' 'TGAT' 'TGCA' 'TGCC' 'TGCG' 'TGCT' 'TGGA' 'TGGC']\n",
      "  ['TTAA' 'TTAC' 'TTAG' 'TTAT' 'TTCA' 'TTCC' 'TTCG' 'TTCT' 'TTGA' 'TTGC']]\n",
      "\n",
      " [['AAAA' 'AAAC' 'AAAG' 'AAAT' 'ACAA' 'ACAC' 'ACAG' 'ACAT' 'AGAA' 'AGAC']\n",
      "  ['AACA' 'AACC' 'AACG' 'AACT' 'ACCA' 'ACCC' 'ACCG' 'ACCT' 'AGCA' 'AGCC']\n",
      "  ['AAGA' 'AAGC' 'AAGG' 'AAGT' 'ACGA' 'ACGC' 'ACGG' 'ACGT' 'AGGA' 'AGGC']\n",
      "  ['AATA' 'AATC' 'AATG' 'AATT' 'ACTA' 'ACTC' 'ACTG' 'ACTT' 'AGTA' 'AGTC']\n",
      "  ['CAAA' 'CAAC' 'CAAG' 'CAAT' 'CCAA' 'CCAC' 'CCAG' 'CCAT' 'CGAA' 'CGAC']\n",
      "  ['CACA' 'CACC' 'CACG' 'CACT' 'CCCA' 'CCCC' 'CCCG' 'CCCT' 'CGCA' 'CGCC']\n",
      "  ['CAGA' 'CAGC' 'CAGG' 'CAGT' 'CCGA' 'CCGC' 'CCGG' 'CCGT' 'CGGA' 'CGGC']\n",
      "  ['CATA' 'CATC' 'CATG' 'CATT' 'CCTA' 'CCTC' 'CCTG' 'CCTT' 'CGTA' 'CGTC']\n",
      "  ['GAAA' 'GAAC' 'GAAG' 'GAAT' 'GCAA' 'GCAC' 'GCAG' 'GCAT' 'GGAA' 'GGAC']\n",
      "  ['GACA' 'GACC' 'GACG' 'GACT' 'GCCA' 'GCCC' 'GCCG' 'GCCT' 'GGCA' 'GGCC']\n",
      "  ['GAGA' 'GAGC' 'GAGG' 'GAGT' 'GCGA' 'GCGC' 'GCGG' 'GCGT' 'GGGA' 'GGGC']\n",
      "  ['GATA' 'GATC' 'GATG' 'GATT' 'GCTA' 'GCTC' 'GCTG' 'GCTT' 'GGTA' 'GGTC']\n",
      "  ['TAAA' 'TAAC' 'TAAG' 'TAAT' 'TCAA' 'TCAC' 'TCAG' 'TCAT' 'TGAA' 'TGAC']\n",
      "  ['TACA' 'TACC' 'TACG' 'TACT' 'TCCA' 'TCCC' 'TCCG' 'TCCT' 'TGCA' 'TGCC']\n",
      "  ['TAGA' 'TAGC' 'TAGG' 'TAGT' 'TCGA' 'TCGC' 'TCGG' 'TCGT' 'TGGA' 'TGGC']\n",
      "  ['TATA' 'TATC' 'TATG' 'TATT' 'TCTA' 'TCTC' 'TCTG' 'TCTT' 'TGTA' 'TGTC']]\n",
      "\n",
      " [['AAAA' 'AACA' 'AAGA' 'AATA' 'ACAA' 'ACCA' 'ACGA' 'ACTA' 'AGAA' 'AGCA']\n",
      "  ['AAAC' 'AACC' 'AAGC' 'AATC' 'ACAC' 'ACCC' 'ACGC' 'ACTC' 'AGAC' 'AGCC']\n",
      "  ['AAAG' 'AACG' 'AAGG' 'AATG' 'ACAG' 'ACCG' 'ACGG' 'ACTG' 'AGAG' 'AGCG']\n",
      "  ['AAAT' 'AACT' 'AAGT' 'AATT' 'ACAT' 'ACCT' 'ACGT' 'ACTT' 'AGAT' 'AGCT']\n",
      "  ['CAAA' 'CACA' 'CAGA' 'CATA' 'CCAA' 'CCCA' 'CCGA' 'CCTA' 'CGAA' 'CGCA']\n",
      "  ['CAAC' 'CACC' 'CAGC' 'CATC' 'CCAC' 'CCCC' 'CCGC' 'CCTC' 'CGAC' 'CGCC']\n",
      "  ['CAAG' 'CACG' 'CAGG' 'CATG' 'CCAG' 'CCCG' 'CCGG' 'CCTG' 'CGAG' 'CGCG']\n",
      "  ['CAAT' 'CACT' 'CAGT' 'CATT' 'CCAT' 'CCCT' 'CCGT' 'CCTT' 'CGAT' 'CGCT']\n",
      "  ['GAAA' 'GACA' 'GAGA' 'GATA' 'GCAA' 'GCCA' 'GCGA' 'GCTA' 'GGAA' 'GGCA']\n",
      "  ['GAAC' 'GACC' 'GAGC' 'GATC' 'GCAC' 'GCCC' 'GCGC' 'GCTC' 'GGAC' 'GGCC']\n",
      "  ['GAAG' 'GACG' 'GAGG' 'GATG' 'GCAG' 'GCCG' 'GCGG' 'GCTG' 'GGAG' 'GGCG']\n",
      "  ['GAAT' 'GACT' 'GAGT' 'GATT' 'GCAT' 'GCCT' 'GCGT' 'GCTT' 'GGAT' 'GGCT']\n",
      "  ['TAAA' 'TACA' 'TAGA' 'TATA' 'TCAA' 'TCCA' 'TCGA' 'TCTA' 'TGAA' 'TGCA']\n",
      "  ['TAAC' 'TACC' 'TAGC' 'TATC' 'TCAC' 'TCCC' 'TCGC' 'TCTC' 'TGAC' 'TGCC']\n",
      "  ['TAAG' 'TACG' 'TAGG' 'TATG' 'TCAG' 'TCCG' 'TCGG' 'TCTG' 'TGAG' 'TGCG']\n",
      "  ['TAAT' 'TACT' 'TAGT' 'TATT' 'TCAT' 'TCCT' 'TCGT' 'TCTT' 'TGAT' 'TGCT']]]\n"
     ]
    }
   ],
   "source": [
    "## What does it look like as character strings\n",
    "mats_strs = np.zeros((3, 16, 16), dtype=\"S4\")\n",
    "mats_strs[0] = pi_arr\n",
    "\n",
    "## fill alternates\n",
    "x = 0\n",
    "for i in xrange(0, 16, 4):\n",
    "    for j in xrange(0, 16, 4):\n",
    "        mats_strs[1, i:i+4, j:j+4] = mats_strs[0, x].reshape(4, 4)\n",
    "        mats_strs[2, i:i+4, j:j+4] = mats_strs[0, x].reshape(4, 4).T\n",
    "        x += 1\n",
    "\n",
    "## print just first 10 columns so that it's more readable\n",
    "print mats_strs[:, :, :10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the matrix rank\n",
    "In this case you can clearly see that the first matrix is the better fit, it has two values that *vanish*, as we expect invariants to do, and thus the rank is 4, where it is 6 for the other two topologies. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mat=0, rank=4\n",
      "mat=1, rank=6\n",
      "mat=2, rank=6\n"
     ]
    }
   ],
   "source": [
    "### Well, first we can just calculate the matrix rank\n",
    "### b/c if they have different rank then we just choose \n",
    "### the lowest one\n",
    "for mat in xrange(3):\n",
    "    print \"mat={}, rank={}\".format(mat, np.linalg.matrix_rank(mats_ints[mat, :, :]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why no SVD? \n",
    "We didn't need it because there were so many empty columns/rows. However, most of the time with real data, and especially large data sets, the matrices will have little or no empty columns, and so the matrices will all be of equal rank. In that case we need to use SVD to calculate which invariants approach most closely to zero from the SVD scores. \n",
    "\n",
    "Here we calculate the squared sum of SVD scores above matrix rank 4, the minimum rank observed in the data set... (not sure why this is, exactly). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mat=0, rank=4, sum=2.08152798678e-15\n",
      "mat=1, rank=6, sum=14.1421356237\n",
      "mat=2, rank=6, sum=14.1421356237\n"
     ]
    }
   ],
   "source": [
    "for idx in xrange(3):\n",
    "    ## if all \n",
    "    np.linalg.matrix_rank(mats_ints[idx])\n",
    "\n",
    "    ## calculate SVD\n",
    "    scores = np.linalg.svd(mats_ints[idx])[1]\n",
    "    \n",
    "    ## sum square of X smallest scores\n",
    "    print \"mat={}, rank={}, sum={}\"\\\n",
    "           .format(idx, \n",
    "                   np.linalg.matrix_rank(mats_ints[idx]),\n",
    "                   np.sqrt(np.sum(scores[4:]**2)))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### There we have it, the first topology (12|34) is best."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}