pso paper implementation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0e07edd6-04c5-4071-b5c8-bc7f050ca5e2",
   "metadata": {},
   "source": [
    "Note to instructors: instructions for running our code are at the bottom of this notebook.\n",
    "All code chunks that were executed during testing of our algorithm were silence with the `ESC + r` hotkey. They can be re-enabled by clicking a chunk and pressing `ESC + y`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58f9f73c-0322-4fea-adfc-5ec783c0cd5e",
   "metadata": {},
   "source": [
    "## Add and load packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ecea89c-4613-4586-b36f-2c060c989670",
   "metadata": {},
   "outputs": [],
   "source": [
    "# install packages\n",
    "import Pkg\n",
    "Pkg.add(\"LinearAlgebra\")\n",
    "Pkg.add(\"Random\")\n",
    "Pkg.add(\"BenchmarkTools\")\n",
    "Pkg.add(\"StaticArrays\")\n",
    "Pkg.add(\"LinearAlgebra\")\n",
    "Pkg.add(\"Distributions\")\n",
    "Pkg.add(\"StatsBase\")\n",
    "Pkg.add(\"Profile\")\n",
    "Pkg.add(\"StatProfilerHTML\")\n",
    "Pkg.add(\"Test\")\n",
    "Pkg.add(\"DataFrames\")\n",
    "Pkg.add(\"Gadfly\")\n",
    "Pkg.add(\"Profile\")\n",
    "Pkg.add(\"PProf\")\n",
    "Pkg.add(\"CSV\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e36f5642",
   "metadata": {},
   "outputs": [],
   "source": [
    "import LinearAlgebra\n",
    "import Random\n",
    "import BenchmarkTools\n",
    "import StaticArrays\n",
    "import LinearAlgebra\n",
    "import Distributions\n",
    "import StatsBase\n",
    "import StatProfilerHTML\n",
    "import Test\n",
    "import CSV\n",
    "using DataFrames\n",
    "using Gadfly\n",
    "using Profile\n",
    "using PProf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71514c58-a48d-47bd-885b-7fcce28de9d7",
   "metadata": {},
   "source": [
    "# Useful Functions for Testing\n",
    "This section contains some basic functionality for generating test data and calculating scores of motifs based on starting positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91d54e65-6664-43b5-b79e-b1602878e80c",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "author: gtb\n",
    "Reads in a fasta file.\n",
    "#### Arguments\n",
    "- `fname`: name of the file to be read\n",
    "\"\"\"\n",
    "# helper function to read in the DNA sequences\n",
    "function ReadInputs(fname)\n",
    "    lines = readlines(fname)\n",
    "    \n",
    "    # collect sequences into here\n",
    "    DNA = Vector{Vector{Char}}()\n",
    "    \n",
    "    # start collecting the first sequence\n",
    "    current = Vector{Char}()\n",
    "    \n",
    "    # go ahead and pop first line\n",
    "    this_line = popat!(lines, 1)\n",
    "    \n",
    "    while length(lines) >= 1\n",
    "        this_line = popat!(lines, 1)\n",
    "        if this_line[1] == '>'\n",
    "            # add the sequence to DNA\n",
    "            push!(DNA, current)\n",
    "            current = Vector{Char}()\n",
    "        else\n",
    "            current = vcat(current, collect(this_line))\n",
    "        end \n",
    "    end\n",
    "    \n",
    "    push!(DNA, current)\n",
    "    # now read the other N lines and collect the strings into a matrix of characters\n",
    "\n",
    "    return(DNA)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5268d27-aea7-495e-a46e-79fe00a8c9da",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "author: gtb\n",
    "Reads in inputs given in the format Dr. Heber posted on Piazza.\n",
    "#### Arguments\n",
    "- `fname`: name of the file to be read\n",
    "\"\"\"\n",
    "# helper function to read in the DNA sequences\n",
    "function ReadInputs_class(fname)\n",
    "    \n",
    "    lines = readlines(fname)\n",
    "    \n",
    "    # collect sequences into here\n",
    "    DNA = Vector{Vector{Char}}()\n",
    "    \n",
    "    # read in the parameters from the first line\n",
    "    lines1 = parse.(Int64, split(popat!(lines, 1)))\n",
    "    \n",
    "    k = lines1[1]\n",
    "    t = lines1[2]    \n",
    "    \n",
    "    for line in lines\n",
    "        push!(DNA, collect(line))\n",
    "    end\n",
    "    # now read the other N lines and collect the strings into a matrix of characters\n",
    "\n",
    "    return(k, DNA)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5421e7dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb\n",
    "Generate a random sequence of length `Length`. Returns a vector of type `Char`.\n",
    "...\n",
    "#### Arguments\n",
    "- `Length`: the length of the sequence to generate using the standard ACGT alphabet\n",
    "...\n",
    "\"\"\"\n",
    "function GenerateSequence(Length)\n",
    "    return collect(Random.randstring(\"ACGT\", Length))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e6553e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb\n",
    "Generate mutliple sequences, each of a given length. Calls `GenerateSequence`. Return a vector of vectors of type `Char`.\n",
    "...\n",
    "#### Arguments\n",
    "- `NumberOfSequences`: the number of sequences to generate\n",
    "- `Length`: the number of nucleotides in each DNA sequences to be generated\n",
    "...\n",
    "\"\"\"\n",
    "function GenerateSequences(NumberOfSequences, Length)\n",
    "    return Sequences = map(i -> GenerateSequence(Length), 1:NumberOfSequences)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ca9ed0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb\n",
    "Mutate a DNA sequence to have a certain distance from the input sequence. Expects and input and output of type `Vec{Char}`.\n",
    "...\n",
    "#### Arguments\n",
    "- `Sequence`: the sequence that will be mutated.\n",
    "- `Distance`: the number of nucleotides in the sequence to change.\n",
    "...\n",
    "\"\"\"\n",
    "function Mutate!(Sequence, Distance) # distance is no. of mutations\n",
    "    # return the sequence if the Distance is 0\n",
    "    if Distance == 0\n",
    "        return Sequence\n",
    "    end\n",
    "    # choose the sites to mutate\n",
    "    posToMutate = StatsBase.sample(1:length(Sequence), Distance, replace = false)\n",
    "    \n",
    "    # pick the letters that can be used at each position\n",
    "    basesAtPositions = [string(i) for i in Sequence[posToMutate]]\n",
    "\n",
    "    lettersToFill = [\"ACGT\" for i in 1:length(Sequence)]\n",
    "    \n",
    "    # figure out what letters each positions can be changed to\n",
    "    lettersToFill = map((x, y) -> replace(x, y => \"\"), lettersToFill, basesAtPositions)\n",
    "    \n",
    "    # select one character from each of these positions and place them where they should be\n",
    "    # use only to cast the string to a character\n",
    "    # https://stackoverflow.com/questions/59946081/julia-convert-string-to-char-or-convert-arraysubstringstring-1-to-char\n",
    "    Sequence[posToMutate] = map(x -> only(Random.randstring(x, 1)), lettersToFill)\n",
    "    \n",
    "    return Sequence\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc4c424f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb\n",
    "Generate test data for testing the (l,d) planted motif problem. Returns the:\n",
    "- `motif` as type `Vec{Char}`, \n",
    "- starting positions of the planted motifs in each sequence `motif_starts` as `Vec{Int64}`,\n",
    "- the actual planted motif including mutations as `motifs_implanted` as `Vec{Vec{Char}}`,\n",
    "- the `sequences` of nucleotides including the motifs as `Vec{Vec{Char}}`.\n",
    "...\n",
    "#### Arguments\n",
    "- `NumberOfSequences`: the number of DNA sequences to produce\n",
    "- `LengthMotif`: the length of the motif to plant\n",
    "- `LengthSequences`: the length of each of the DNA sequences the motif will be planted into\n",
    "- `Distance`: the hamming distance of each planted sequences from the consensus motif\n",
    "...\n",
    "\"\"\"\n",
    "function GenerateTestData_ld(NumberOfSequences, LengthMotif, LengthSequences, Distance)\n",
    "    # do some error checking to make sure the values provided are valid\n",
    "    # specifically, Disance <= LengthMotif\n",
    "    if Distance > LengthMotif\n",
    "        error(\"The Distance if larger than the LengthMotif.\")\n",
    "    end\n",
    "\n",
    "    # LengthMotif <= LengthSequences\n",
    "    if LengthMotif > LengthSequences\n",
    "        error(\"The LengthMotif is longer than the LengthSequences\")\n",
    "    end\n",
    "    \n",
    "    # NumberOfSequences >= 2\n",
    "    if NumberOfSequences <= 1\n",
    "        error(\"The NumberOfSequences is 1, which is too small for motif detection\")\n",
    "    end\n",
    "\n",
    "    # generate the input sequences\n",
    "    sequences = GenerateSequences(NumberOfSequences, LengthSequences)\n",
    "    \n",
    "    # make the motif\n",
    "    motif = GenerateSequence(LengthMotif)\n",
    "    \n",
    "    # generate a mutated motif for each sequence to implant\n",
    "    motifs = [copy(motif) for i in 1:NumberOfSequences]\n",
    "    \n",
    "    # broadcast the mutated motif\n",
    "    motifs .= Mutate!.(motifs, Distance)\n",
    "    \n",
    "    # now implant the mutated motifs and record where we implant it\n",
    "    motifStarts = rand(1:(LengthSequences - LengthMotif + 1), NumberOfSequences)\n",
    "    \n",
    "    # now place the motifs into the sequences vector\n",
    "    for i in 1:length(sequences)\n",
    "        sequences[i][motifStarts[i]:(motifStarts[i] + LengthMotif - 1)] = motifs[i]\n",
    "    end\n",
    "    \n",
    "    # now return all the things we may want later\n",
    "    return (motif = motif, motifs_starts = motifStarts, motifs_implanted = motifs, sequences = sequences)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b262b0f-9508-46a7-b94a-22db8e197ca6",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb\n",
    "Generate test data for testing the (l,d) planted motif problem. Returns the:\n",
    "- `motif` as type `Vec{Char}`, \n",
    "- starting positions of the planted motifs in each sequence `motif_starts` as `Vec{Int64}`,\n",
    "- the actual planted motif including mutations as `motifs_implanted` as `Vec{Vec{Char}}`,\n",
    "- the `sequences` of nucleotides including the motifs as `Vec{Vec{Char}}`.\n",
    "...\n",
    "#### Arguments\n",
    "- `NumberOfSequences`: the number of DNA sequences to produce\n",
    "- `LengthMotif`: the length of the motif to plant\n",
    "- `LengthSequences`: the length of each of the DNA sequences the motif will be planted into\n",
    "- `Distance`: the hamming distance of each planted sequences from the consensus motif\n",
    "- `MissingPerc\n",
    "...\n",
    "\"\"\"\n",
    "function GenerateTestData_ld_somemissing(NumberOfSequences, LengthMotif, LengthSequences, Distance, MissingPerc)\n",
    "    # do some error checking to make sure the values provided are valid\n",
    "    # specifically, Disance <= LengthMotif\n",
    "    if Distance > LengthMotif\n",
    "        error(\"The Distance if larger than the LengthMotif.\")\n",
    "    end\n",
    "\n",
    "    # LengthMotif <= LengthSequences\n",
    "    if LengthMotif > LengthSequences\n",
    "        error(\"The LengthMotif is longer than the LengthSequences\")\n",
    "    end\n",
    "    \n",
    "    # NumberOfSequences >= 2\n",
    "    if NumberOfSequences <= 1\n",
    "        error(\"The NumberOfSequences is 1, which is too small for motif detection\")\n",
    "    end\n",
    "\n",
    "    # generate the input sequences\n",
    "    sequences = GenerateSequences(NumberOfSequences, LengthSequences)\n",
    "    \n",
    "    # make the motif\n",
    "    motif = GenerateSequence(LengthMotif)\n",
    "    \n",
    "    # generate a mutated motif for each sequence to implant\n",
    "    motifs = [copy(motif) for i in 1:NumberOfSequences]\n",
    "    \n",
    "    # broadcast the mutated motif\n",
    "    motifs .= Mutate!.(motifs, Distance)\n",
    "    \n",
    "    # now implant the mutated motifs and record where we implant it\n",
    "    motifStarts = rand(1:(LengthSequences - LengthMotif + 1), NumberOfSequences)\n",
    "    \n",
    "    # figure out which sequences to insert into\n",
    "    insert_list = StatsBase.sample(1:length(sequences), floor(Int64, (1 - MissingPerc) * NumberOfSequences), replace = false)\n",
    "    # now place the motifs into the sequences vector\n",
    "    for i in insert_list\n",
    "        sequences[i][motifStarts[i]:(motifStarts[i] + LengthMotif - 1)] = motifs[i]\n",
    "    end\n",
    "    \n",
    "    # now return all the things we may want later\n",
    "    return (motif = motif, motifs_starts = motifStarts, motifs_implanted = motifs, sequences = sequences)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0013b4af",
   "metadata": {},
   "source": [
    "# PSO Helper Functions\n",
    "Helper functions that are useful for the implementation described in our template paper."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d83e4afc-3347-49f5-b577-5e3802638e83",
   "metadata": {},
   "source": [
    "## Basic Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30f51d1e-041b-40e7-b3a4-d5eee0f660dd",
   "metadata": {},
   "source": [
    "### make_profile\n",
    "This function is used to turn a set of sequences into a sequence profile. Pseudo_counts are enabled by default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94da59a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: ks; mod gtb 11/10  \n",
    "\n",
    "Returns the profile probabilities for the given inputs. Normalized frequencies are given in the order A, C, G, T.\n",
    "\n",
    "#### Arguments\n",
    "- `k`: length of motifs\n",
    "- `t`: number of sequences\n",
    "- `arr`: sequences to operate over, given as a `Vec{Vec{Char}}`\n",
    "- `pseudo_counts`: optional argument for whether or not pseudo counts are used. Defaults to `true`.\n",
    "\"\"\"\n",
    "function make_profile(k, t, arr; pseudo_counts=true) \n",
    "    res = zeros(4,k)\n",
    "    if pseudo_counts==false\n",
    "        for i in 1:k\n",
    "            c = Dict{Char, Int64}('A'=>0,'C'=>0,'G'=>0,'T'=>0)\n",
    "            for j in 1:t\n",
    "                c[arr[j][i]]+=1\n",
    "            end\n",
    "            res[1,i] = c['A']/t\n",
    "            res[2,i] = c['C']/t\n",
    "            res[3,i] = c['G']/t\n",
    "            res[4,i] = c['T']/t\n",
    "\n",
    "        end\n",
    "    else\n",
    "        for i in 1:k\n",
    "            c = Dict{Char, Int64}('A'=>1,'C'=>1,'G'=>1,'T'=>1)\n",
    "            for j in 1:t\n",
    "                c[arr[j][i]]+=1\n",
    "            end\n",
    "            res[1,i] = (c['A'])/(4+t)\n",
    "            res[2,i] = (c['C'])/(4+t)\n",
    "            res[3,i] = (c['G'])/(4+t)\n",
    "            res[4,i] = (c['T'])/(4+t)\n",
    "\n",
    "        end\n",
    "    end\n",
    "    return res\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e902818c-f503-46bd-b227-69ea74ef3b13",
   "metadata": {},
   "source": [
    "### background_frequency\n",
    "This function is used to find the background frequencies for each nucleotide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9729e484-104f-4707-91eb-ebb5765f0958",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: ks  \n",
    "\n",
    "Returns the background frequency of each nucleotide as a dictionary.\n",
    "\n",
    "#### Arguments\n",
    "- `sequence`: the sequence as `Vec{Char}` that will be counted.\n",
    "- `prop`: optional parameter - will convert values to probabilities\n",
    "\"\"\"\n",
    "function background_frequency(sequence; prop=false)\n",
    "    bg = StatsBase.countmap(sequence)\n",
    "    # add keys if they are not present\n",
    "    if !('A' in keys(bg))\n",
    "        bg['A'] = 0\n",
    "    end\n",
    "    if !('C' in keys(bg))\n",
    "        bg['C'] = 0\n",
    "    end\n",
    "    if !('G' in keys(bg))\n",
    "        bg['G'] = 0\n",
    "    end\n",
    "    if !('T' in keys(bg))\n",
    "        bg['T'] = 0\n",
    "    end\n",
    "    if prop\n",
    "        bg_prop = Dict{Char, Float64}()\n",
    "        for k in keys(bg)\n",
    "            bg_prop[k] = bg[k]/length(sequence)\n",
    "        end\n",
    "        return(bg_prop)\n",
    "    else\n",
    "        return(bg)\n",
    "    end\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0bbc8ba-e548-49dd-a6ec-01c4cf288dc1",
   "metadata": {},
   "source": [
    "### get_index\n",
    "This function is used to map the nucleotides to the integers. Useful to have at compile time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c807455d-e060-4c82-b7e0-5e5bfc3ff304",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: ks, mod gtb 11/10  \n",
    "\n",
    "Helper function that maps the nucleotides `[A, C, G, T]` to the integers `1:4`. Useful for accessing certain arrays.\n",
    "\n",
    "TODO: check a valid nucleotide is being passed\n",
    "\n",
    "#### Arguments\n",
    "- `c`: the nucleotide character that is being accessed.\n",
    "\"\"\"\n",
    "function get_index(c)\n",
    "    #lobal d = Dict('A'=>1,'C'=>2,'G'=>3,'T'=>4)\n",
    "    if c == 'A'\n",
    "        return 1\n",
    "    elseif c == 'C'\n",
    "        return 2\n",
    "    elseif c == 'G'\n",
    "        return 3\n",
    "    elseif c == 'T'\n",
    "        return 4\n",
    "    else\n",
    "        print(\"ERROR!!\")\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab3b3fc2-1ef2-4a93-a1d4-63878c4ca4e8",
   "metadata": {},
   "source": [
    "### consensus_score\n",
    "This function uses the position weight matrix to calculate a score for a given set of starting positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4080a10-3cae-4b88-a370-de695a8b30c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: ks, mod gtb 11/11\n",
    "\n",
    "Helper function to calculate the consensus score given a set of motifs and the background frequencies of the nucleotides.\n",
    "\n",
    "#### Arguments \n",
    "- `matches`: obtained motif-matching substrings for all sequences\n",
    "- `background_frequency`: probability of each base in the background\n",
    "\"\"\"\n",
    "function consensus_score(matches, background_frequency)\n",
    "    # build the map from numbers to nucleotides\n",
    "    int_char = Dict{Int8, Char}(1=>'A',2=>'C',3=>'G',4=>'T')\n",
    "    \n",
    "    # determine the length of the motif\n",
    "    k = length(matches[1])\n",
    "    t = length(matches)\n",
    "    \n",
    "    # build a profile for each column\n",
    "    profile = make_profile(k, t, matches; pseudo_counts=false)\n",
    "    \n",
    "    # in case of ties it sticks to A > C > G > T\n",
    "    # generate a vector for storing the most common character in each column\n",
    "    \n",
    "    # find the most common value in each column\n",
    "    # cmotif = StaticArrays.SVector{k, Char}([int_char[i] for i in vec(mapslices(i -> findmax(i)[2], profile; dims = 1))])\n",
    "    cmotif = StaticArrays.SVector{k, Char}([int_char[i] for i in [ch[1] for ch in findmax(profile, dims = 1)[2]]])\n",
    "    \n",
    "    # now sum over the calculation in the paper\n",
    "    s = 0\n",
    "    for i in 1:k\n",
    "        for j in 1:4\n",
    "            idx_char = int_char[j]\n",
    "            toadd = profile[j, i] * log(2, profile[j, i]/background_frequency[idx_char])\n",
    "            # add nothing if we get NaN, since this corresponds to no matches\n",
    "            if !isnan(toadd)\n",
    "                s += toadd\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return s, cmotif\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d5a7ace-21fd-4635-b9bc-99501c649bcd",
   "metadata": {},
   "source": [
    "### weights\n",
    "This function is used to turn nucleotide frequencies into weights, similar to how our template paper describes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "671864c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: ks  \n",
    "\n",
    "Calculates the weights for each nucleotide based on their frequencies. Returns these values as a `Dict{Char -> Float64}`.\n",
    "\n",
    "#### Arguments\n",
    "- `background_frequency`: the background frequencies of each nucleotide from `background_frequency`\n",
    "\"\"\"\n",
    "function weights(background_frequency)\n",
    "    ma = max(1/background_frequency['A'],1/background_frequency['C'],\n",
    "        1/background_frequency['G'],1/background_frequency['T'])\n",
    "    \n",
    "    d = zeros(Float64, 4)\n",
    "    \n",
    "    d[1] = (1/background_frequency['A'])/ma\n",
    "    d[2] = (1/background_frequency['C'])/ma\n",
    "    d[3] = (1/background_frequency['G'])/ma\n",
    "    d[4] = (1/background_frequency['T'])/ma\n",
    "    \n",
    "    # failsafe: if any are NaN, return all ones\n",
    "    if any(isnan.(d))\n",
    "        return ones(Float64, 4)\n",
    "    else\n",
    "        return d\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "954b88d5-3f51-4d71-a99e-0c069e51c74d",
   "metadata": {},
   "source": [
    "### update_motif\n",
    "Given the current particle, its personal best, and the global best, this function proposes a new particle that will then be used for search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eaff733-07b4-4997-82cf-d609088e5f7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: ks, mod gtb 11/11  \n",
    "\n",
    "Update the current nucleotide for a particle given the personal best, global best, and the weighting vectors for the vectors, `c`, and nucleotides, `w`.\n",
    "\n",
    "#### Arguments\n",
    "- `current_seq`: the current motif sequence for the particle\n",
    "- `pbest_seq`: the motif sequence that corresponds to the best this particle has ever been\n",
    "- `gbest_seq`: the motif sequence corresponding to the best any particle has achieved\n",
    "- `weights`: a vector of weights for each nucleotide, based on their frequencies\n",
    "- `background_frequency`: probability of each base in the background\n",
    "\"\"\"\n",
    "function update_motif(current_seq, pbest_seq, gbest_seq, weights, background_frequency; deterministic=false)\n",
    "    k = length(current_seq)\n",
    "    \n",
    "    if deterministic\n",
    "        Random.seed!(101)\n",
    "        scale_vals = ones(k, 4) .* permutedims(repeat(weights', k))\n",
    "        rand_seq = collect(Random.randstring(['A','C','G','T'], k))\n",
    "    else\n",
    "        rand_seq = collect(Random.randstring(['A','C','G','T'], k))\n",
    "        scale_vals = rand(4, k) .* permutedims(repeat(weights', k))\n",
    "    end\n",
    "\n",
    "    # bind all of the sequences together\n",
    "    seq_bind = permutedims(reduce(hcat, [current_seq, pbest_seq, gbest_seq, rand_seq]))\n",
    "    \n",
    "    # now get the cumulative sum in each column so we can sample a probability to figure out which base to keep\n",
    "    # now return the \"best\" nucleotide for that position based on our rules\n",
    "    return(seq_bind[vec(findmax(scale_vals; dims = 1)[2])])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e916aa0-49a8-4777-9934-60fc8bd3fde6",
   "metadata": {},
   "source": [
    "# Our Creative Contributions\n",
    "This section contains most of our group's creative contributions. We use pre-computation, indexing, and optimistic score-thresholding to reduce the number of wasteful computations. How these functions work may be more obvious by looking at the test cases."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc75f477-e02b-48b5-be40-550d245b5157",
   "metadata": {},
   "source": [
    "## Precomputation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df9c3797-c0d4-4dcb-8317-f98d41fc70e9",
   "metadata": {},
   "source": [
    "### precompute_matchscores\n",
    "This function is used to avoid having to do expensive log calculations thousands of times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edfcc645-23c5-4d0e-8b1c-90c838eb0081",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Precompute the log match scores. Returns a dictionary whose keys are (MOTIF LETTER, SEQUENCE LETTER)\n",
    "\n",
    "#### Arguments\n",
    "- `background_frequency`: the background frequency for each nucleotide\n",
    "\"\"\"\n",
    "function precompute_matchscores(background_frequency)\n",
    "    pc_dict = Dict{Tuple{Char, Char}, Float64}()\n",
    "    for i in ['A', 'C', 'G', 'T']\n",
    "        for j in ['A', 'C', 'G', 'T']\n",
    "            if i == j\n",
    "                pc_dict[(i, j)] = 1 + log(4, (0.25/background_frequency[i]))\n",
    "            else\n",
    "                pc_dict[(i, j)] = log(4, 0.25/sqrt(background_frequency[i] * background_frequency[j]))\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return(pc_dict)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4db03d6-793c-4432-8f1a-8c49fc6d2d47",
   "metadata": {},
   "source": [
    "### precompute_hash\n",
    "This function creates an indexible dictionary for trying to find the substring most similar to a proposed motif."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5423a359-d1c5-4180-9540-fdcb397a88dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb 11/23\n",
    "\n",
    "Returns a dictionary with the following structure:\n",
    "    Level 1: positions 1:k\n",
    "    Level 2: characters A, C, G, T\n",
    "    Level 3: starting positions with this value\n",
    "\n",
    "#### Arguments\n",
    "- `seq`: the sequence to be hashed\n",
    "- `k`: the length of the pattern to be hashed\n",
    "\"\"\"\n",
    "function precompute_hash(seq, k)\n",
    "    # make a dictionary with\n",
    "    out_dict = Dict{Int64, Dict{Char, BitVector}}()\n",
    "    \n",
    "    # get all chars\n",
    "    allchars = permutedims(reduce(hcat, [seq[i:(i+k-1)] for i in 1:(length(seq)-k+1)]))\n",
    "    \n",
    "    for pos in 1:k\n",
    "        intermediate_dict = Dict{Char, BitVector}()\n",
    "        for ch in ['A', 'C', 'G', 'T']\n",
    "            intermediate_dict[ch] = allchars[:,pos] .== ch\n",
    "        end\n",
    "        out_dict[pos] = intermediate_dict\n",
    "    end\n",
    "    return(out_dict)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b29a6c4-3439-4f74-b0e0-e48b9b102771",
   "metadata": {},
   "source": [
    "### precompute_motifscores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14a64f22-b53f-4072-80b5-f0e002f9624f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author gtb 11/26\n",
    "\n",
    "Precomputes the \"phases\" of the motif, so we know what the fitness will be along the entire sequence, just by knowing which base is in which position.\n",
    "\n",
    "Will store these values in a Dictionary with the followign structure:\n",
    "    {Motif Pos -> {{Char in Motif} -> {Vector of Float}}\n",
    "\n",
    "#### Arguments\n",
    "- `seq`: the sequence to be hashed\n",
    "- `seq_hashed`: the sequence that has already been hashed, with the BitVector method\n",
    "- `k` the length of the motif\n",
    "- `matchscores` precomputed values from precompute_matchscores\n",
    "\"\"\"\n",
    "function precompute_motifscores(seq, seq_hashed, k, matchscores)\n",
    "    ms_dict = Dict{Int64, Dict{Char, Vector{Float64}}}()\n",
    "    # iterate over the positions in the motif\n",
    "    for m in 1:k\n",
    "        # create a dict to store value for this match position\n",
    "        intermediate_dict = Dict{Char, Vector{Float64}}()\n",
    "        # add a key for if the base in the motif is A, C, G, or T\n",
    "        for ch1 in ['A', 'C', 'G', 'T']\n",
    "            # set it to all zeroes to start\n",
    "            intermediate_dict[ch1] = zeros(Float64, length(seq) - k + 1)\n",
    "            # now add in the corresponding score for each start position in the sequence\n",
    "            for ch2 in ['A', 'C', 'G', 'T']\n",
    "                mscore = matchscores[(ch1, ch2)]\n",
    "                intermediate_dict[ch1][seq_hashed[m][ch2]] .+= mscore\n",
    "            end\n",
    "        end\n",
    "        # add the intermediate dict to the big dict\n",
    "        ms_dict[m] = intermediate_dict\n",
    "    end\n",
    "    return(ms_dict)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5c2ea0e-95fe-409b-8963-10c78d788c8c",
   "metadata": {},
   "source": [
    "## Faster Searches"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "156e7e87-f49e-432a-8303-6365b8a8f067",
   "metadata": {},
   "source": [
    "### best_match_efficient\n",
    "This function finds the substring that matches most closely to the proposed motif. It only does addition and indexing operations, making it have very low overhead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2db9987-24b1-4a63-917a-de329a478725",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb 11/23\n",
    "\n",
    "Returns the score of the closest matching subsequence in `seq` to `motif`, in a more efficient way.\n",
    "\n",
    "#### Arguments\n",
    "- `seq`: the sequence to search\n",
    "- `motifscore_hashed`: the hashed scores from precompute_motifscores\n",
    "- `motif`: the motif to search for\n",
    "- `matchschores`: precomputed values from precompute_matchscores\n",
    "- `scores`: the matrix scores will be stored in\n",
    "\"\"\"\n",
    "function best_match_efficient(seq, motifscore_hashed, motif, matchscores, scores)  \n",
    "    # iterate over the positions in the motif\n",
    "    for m in 1:length(motif)\n",
    "        @views scores .+= motifscore_hashed[m][motif[m]]\n",
    "    end\n",
    "    \n",
    "    return(findmax(scores))        \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b242f805-9f8f-466f-a790-932ea466f508",
   "metadata": {},
   "source": [
    "### best_possible\n",
    "This function figures out what the highest possible score you could get given the background frequencies, current profile, the number of sequences examined, and the total number of sequences there are to examine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88b09f60-5143-48de-af04-6c672926da70",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb\n",
    "\n",
    "Finds the maximum fitness that could be achieved with a given motif profile\n",
    "\n",
    "#### Arguments\n",
    "- `probs_matrix`: the background nucleotide frequencies, in a special matrix format\n",
    "- `profile`: the profile up until this point\n",
    "- `nseq_sofar`: the number of sequences scored so far\n",
    "- `nseq`: the number of sequences being scored total\n",
    "- `matrix_forfindmax`: the matix that will be used for find max (4 x length motif)\n",
    "- `bestpossible_fitness`: a pre-allocated vector Float64 that is the length of the motif\n",
    "\"\"\"\n",
    "function best_possible(probs_matrix, profile, nseq_sofar, nseq, matrix_forfindmax, bestpossible_fitness)\n",
    "    # iterate over the columns of the profile\n",
    "    for j in 1:size(profile, 2)\n",
    "        # place the current values into matrix_forfindmax\n",
    "        @inbounds matrix_forfindmax .= profile[1:4, j]\n",
    "        \n",
    "        # now add in the values along nseq\n",
    "        matrix_forfindmax[LinearAlgebra.diagind(matrix_forfindmax)] .+= nseq - nseq_sofar\n",
    "        \n",
    "        # convert to probabilities\n",
    "        matrix_forfindmax ./= nseq\n",
    "        \n",
    "        # now do the actual math\n",
    "        matrix_forfindmax .= matrix_forfindmax .* log2.(matrix_forfindmax ./ probs_matrix)\n",
    "        \n",
    "        matrix_forfindmax[isnan.(matrix_forfindmax)] .= 0\n",
    "\n",
    "        # sum over the columns\n",
    "        bestpossible_fitness[j] = maximum(sum(matrix_forfindmax; dims = 1))\n",
    "    end\n",
    "    \n",
    "    return(sum(bestpossible_fitness))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd045e9d-da35-477c-8d89-ac49dec1d654",
   "metadata": {},
   "source": [
    "### optimistic_search\n",
    "This function will look for a given motif in each of the provided sequences. It's main advantage is that it quits after it realizes that it cannot beat the current best score for the particle. It does this calculation with the help of `best_possible`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a18f54a-b88f-4a35-bf2b-ef41606400a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: gtb\n",
    "\n",
    "Searches using a motif consensus sequence until the optimistic score cannot reach pbest\n",
    "\n",
    "#### Arguments\n",
    "- `seqs`: the sequences that are to be searched\n",
    "- `motif`: the motif sequence to search for\n",
    "- `background_frequency`: the background nucleotide frequencies\n",
    "- `pbest_thresh`: the current best fitness of the particle\n",
    "- `motifscore_hashed`: the hashed scores, a dict that makes it faster to search\n",
    "- `precomp`: the precomputed scoring calculations (no more logs!)\n",
    "\"\"\"\n",
    "function optimistic_search(seqs, motif, background_frequency, pbest_thresh, motifscore_hashed, precomp; vanilla = false)\n",
    "    best_match_starts = zeros(Int64, length(seqs))\n",
    "\n",
    "    # start a profile we will add on to\n",
    "    search_profile = zeros(Float64, 4, length(motif))\n",
    "    \n",
    "    optimal = zeros(Float64, 1)\n",
    "    \n",
    "    # convert the background frequencies into a matrix\n",
    "    probs_matrix = permutedims(repeat([background_frequency['A'] background_frequency['C'] background_frequency['G'] background_frequency['T']], 4))\n",
    "    \n",
    "    # create a 4x4 matrix we will use for each column; we will keep putting values in here\n",
    "    mfm = zeros(Float64, 4, 4)\n",
    "    \n",
    "    # create a vector to store the current fitness in each column of the profile matrix\n",
    "    bpf = zeros(Float64, size(search_profile, 2))\n",
    "    \n",
    "     # create a l - k + 1 vector for storing the scores\n",
    "    scores = zeros(Float64, length(seqs[1]) - length(motif) + 1)\n",
    "    \n",
    "    for seq in 1:length(seqs)\n",
    "        # make sure scores is set to zero before passing to best_match_efficient\n",
    "        scores .= 0\n",
    "        \n",
    "        # find the best match subsequence\n",
    "        best_match_starts[seq] = best_match_efficient(seqs[seq], motifscore_hashed[seq], motif, precomp, scores)[2]\n",
    "        \n",
    "        # iterate along the best match string\n",
    "        for ch in 1:length(motif)\n",
    "            search_profile[get_index(seqs[seq][best_match_starts[seq] + ch - 1]), ch] += 1 \n",
    "        end\n",
    "        \n",
    "        # skip the break condition if we are doing vanilla search\n",
    "        if vanilla\n",
    "            continue\n",
    "        end\n",
    "        \n",
    "        optimal[1] = best_possible(probs_matrix, search_profile, seq, length(seqs), mfm, bpf)\n",
    "\n",
    "        # find the optimistic score\n",
    "        if optimal[1] < pbest_thresh\n",
    "            # return some condition that says this particle isn't going to beat the personal best\n",
    "            return(best_match_starts)\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    # return at the end\n",
    "    return(best_match_starts)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "140c913b-cf21-4e64-bdb4-4d37130141a3",
   "metadata": {},
   "source": [
    "# PSO Implementation\n",
    "This program actually runs the particle swarm optimization.\n",
    "It has complexity approximately equal to:  \n",
    "`# of initializations` $\\times$ `# of particles` $\\times$ `# of iterations` $\\times$ `(sequence lengths - motif length + 1)` $\\times$ `# of sequences - 1`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15eef2ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: ks  \n",
    "\n",
    "PSO Motif is the main function for this project. It takes some set of parameters, and then calls other functions to do the actual work.\n",
    "\n",
    "#### Arguments\n",
    "- `sequences`: the input DNA sequences, stored as a `Vec{Vec{Char}}`\n",
    "- `motiflen`: the motif length to look for\n",
    "- `max_reset`: the number of times to re-intialize\n",
    "- `max_iteration`: the number of times to update each particle\n",
    "- `nparticles`: the number of particles to track\n",
    "\"\"\"\n",
    "function pso_motif(sequences, motiflen, max_reset, max_iteration, nparticles; includestarts = false)\n",
    "    ## hash all the sequences so they are easier to search\n",
    "    sh = precompute_hash.(sequences, (motiflen,))\n",
    "    \n",
    "    # set the particles fitness to negative infinity to start\n",
    "    final_fitness = -Inf16\n",
    "    \n",
    "    # create a vector to store the final motif in\n",
    "    final_motif = []\n",
    "    \n",
    "    # and the starting positions\n",
    "    final_starts = []\n",
    "    \n",
    "    # set the weight parameters\n",
    "    # current, pbest, gbest, random\n",
    "    c = StaticArrays.SVector{4, Float64}([1.0, 1.0, 1.0, 1.0])\n",
    "    \n",
    "    # calculate the background frequencies\n",
    "    bg = background_frequency(reduce(hcat, sequences); prop = true)\n",
    "    \n",
    "    # precompute all the scores for each type of match\n",
    "    ms = precompute_matchscores(bg)\n",
    "    \n",
    "    # hash the motif scores as well\n",
    "    msh = precompute_motifscores.(sequences, sh, (motiflen,), (ms,))\n",
    "    \n",
    "    # and their weights\n",
    "    w = weights(bg)\n",
    "    \n",
    "    # initialize some arrays!\n",
    "    particle_startingpos = zeros(Int64, nparticles)\n",
    "    particle_startingseq = zeros(Int64, nparticles)\n",
    "    \n",
    "    particles = [StaticArrays.SVector{motiflen, Char}(['A' for i in 1:motiflen]) for i in 1:nparticles]\n",
    "\n",
    "    match_sequences = [['A' for i in 1:motiflen] for i in 1:length(sequences)]\n",
    "    \n",
    "    # pbest will be an array of arrays, corresponding to the best consensus motif each particle has personally seen\n",
    "    pbest = [StaticArrays.SVector{motiflen, Char}(['A' for i in 1:motiflen]) for i in 1:nparticles]\n",
    "    \n",
    "    for i in 1:max_reset \n",
    "        # sample particle starting positions\n",
    "        # and the sequences they come from\n",
    "        \n",
    "        ## TODO: maybe ensure the starting positions result in non-overlaping particles? in terms of consensus motifs\n",
    "        particle_startingpos .= StatsBase.sample(1:(length(sequences[1]) - motiflen + 1), nparticles, replace = true)\n",
    "        particle_startingseq .= StatsBase.sample(1:length(sequences), nparticles, replace = true)\n",
    "\n",
    "        # extract these sequences from the input sequences\n",
    "        for p in 1:nparticles\n",
    "            # add these particles into a static array\n",
    "            particles[p] = StaticArrays.SVector{motiflen, Char}(sequences[particle_startingseq[p]][particle_startingpos[p]:(particle_startingpos[p]+motiflen-1)])\n",
    "            pbest[p] = StaticArrays.SVector{motiflen, Char}(sequences[particle_startingseq[p]][particle_startingpos[p]:(particle_startingpos[p]+motiflen-1)])\n",
    "        end\n",
    "        \n",
    "        # set the personal best fitness for each particle\n",
    "        # and for global to be -Inf\n",
    "        # and create a vector to store the current scores inside of\n",
    "        fitness_current = [-Inf16 for i in 1:nparticles]\n",
    "        fitness_pbest = [-Inf16 for i in 1:nparticles]\n",
    "        fitness_gbest = [-Inf16]\n",
    "        \n",
    "        gbest = [StaticArrays.SVector{motiflen, Char}(['A' for i in 1:motiflen])]\n",
    "            \n",
    "        # iterate until convergance\n",
    "        j = 1\n",
    "        no_updates = 0\n",
    "\n",
    "        match_positions = zeros(Int64, length(sequences))\n",
    "        best_match_saves_starts = [zeros(Int64, length(sequences)) for i in 1:nparticles]\n",
    "        while (j <= max_iteration) & (no_updates <= 5)           \n",
    "            for k in 1:nparticles\n",
    "                #=display(fitness_current)\n",
    "                display(fitness_pbest)\n",
    "                display(fitness_gbest)\n",
    "                display(particles)\n",
    "                display(pbest)\n",
    "                display(gbest[1])\n",
    "                # first, query each sequence to see if we can find the best match to the motif\n",
    "                # do this until we find that our optimistic score cannot be as good as fitness_pbest (and by extension, gbest)\n",
    "                # basically just keep collecting sequences until we're sure we can't do as well as we've done in the past\n",
    "                # then, we can quit trying to update this particle\n",
    "                # basically, this will be a short-circuiting operation=#\n",
    "                match_positions .= optimistic_search(sequences, particles[k], bg, fitness_pbest[k], msh, ms; vanilla = true)\n",
    "\n",
    "                # skip the rest if we aborted the search\n",
    "                if match_positions[length(sequences)] == 0\n",
    "                    continue\n",
    "                end\n",
    "                \n",
    "                for ms in 1:length(sequences)\n",
    "                    #match_sequences[ms] = StaticArrays.SVector{motiflen, Char}(sequences[ms][match_positions[ms]:(match_positions[ms] + motiflen - 1)])\n",
    "                    @views match_sequences[ms] .= sequences[ms][match_positions[ms]:(match_positions[ms] + motiflen - 1)]\n",
    "                end\n",
    "                \n",
    "                # and calculate the consensus score like the paper says - this will become our fitness\n",
    "                fitness_current[k], particles[k] = consensus_score(match_sequences, bg)\n",
    "\n",
    "                if fitness_current[k] > fitness_pbest[k]\n",
    "                    best_match_saves_starts[k] .= copy(match_positions)\n",
    "                    fitness_pbest[k] = copy(fitness_current[k])\n",
    "                    pbest[k] = particles[k]\n",
    "                end\n",
    "                \n",
    "                # update global best if this is the best we've ever had\n",
    "                if fitness_current[k] > fitness_gbest[1]\n",
    "                    final_starts = copy(match_positions)\n",
    "                    fitness_gbest[1] = copy(fitness_current[k])\n",
    "                    gbest[1] = particles[k]\n",
    "                end   \n",
    "            end\n",
    "\n",
    "            # do check shift one in every 10 times\n",
    "            if (j % floor(max_iteration/10)) == 0\n",
    "                for k in 1:nparticles\n",
    "                    starts_toworkfrom = copy(best_match_saves_starts[k])\n",
    "                    for shift in [-5, -4, -3, -2, -1, 1, 2, 3, 4, 5]\n",
    "                        tocheck_matchstarts = copy(starts_toworkfrom)\n",
    "                        tocheck_matchstarts .+= shift\n",
    "                        tocheck_matchstarts[tocheck_matchstarts .< 1] .= 1\n",
    "                        tocheck_matchstarts[tocheck_matchstarts .> (length(sequences[1]) - motiflen + 1)] .= (length(sequences[1]) - motiflen + 1)\n",
    "                        \n",
    "                        for ms in 1:length(sequences)\n",
    "                            @views match_sequences[ms] .= sequences[ms][tocheck_matchstarts[ms]:(tocheck_matchstarts[ms] + motiflen - 1)]\n",
    "                        end\n",
    "                        \n",
    "                        potential_fitness_current, potential_new_part = consensus_score(match_sequences, bg)\n",
    "                        \n",
    "                        if potential_fitness_current > fitness_current[k]\n",
    "                            fitness_current[k] = copy(potential_fitness_current)\n",
    "                            particles[k] = StaticArrays.SVector{motiflen, Char}(potential_new_part)\n",
    "                            best_match_saves_starts[k] = copy(tocheck_matchstarts)\n",
    "                        end\n",
    "                        \n",
    "                        if fitness_current[k] > fitness_pbest[k]\n",
    "                            fitness_pbest[k] = copy(fitness_current[k])\n",
    "                            pbest[k] = particles[k]\n",
    "                        end\n",
    "\n",
    "                        # upate global best if this is the best we've ever had\n",
    "                        if fitness_current[k] > fitness_gbest[1]\n",
    "                            no_updates = 0\n",
    "                            final_starts = copy(tocheck_matchstarts)\n",
    "                            fitness_gbest[1] = copy(fitness_current[k])\n",
    "                            gbest[1] = particles[k]\n",
    "                        end  \n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "            \n",
    "            no_updates += 1\n",
    "            \n",
    "            # broadcast the update operation across all the particles\n",
    "            particles .= update_motif.(particles, pbest, (gbest[1], ), (w, ), (bg, ))\n",
    "\n",
    "            j += 1\n",
    "\n",
    "            # break if >50% of particles are the same\n",
    "            #=\n",
    "            if maximum(values(StatsBase.countmap(pbest))) > floor(nparticles * 0.90)\n",
    "                break\n",
    "            end=#\n",
    "        end\n",
    "        \n",
    "        # wrap up by looking for the motif sequence in all sequences one more time\n",
    "        \n",
    "        #=\n",
    "        display(fitness_current)\n",
    "        display(fitness_pbest)\n",
    "        display(fitness_gbest)\n",
    "        display(particles)\n",
    "        display(pbest)\n",
    "        display(gbest[1])=#\n",
    "        \n",
    "        # get the _actual_ best score using all the data\n",
    "        match_positions = optimistic_search(sequences, gbest[1], bg, -Inf, msh, ms; vanilla = true)\n",
    "\n",
    "        for ms in 1:length(sequences)\n",
    "            @views match_sequences[ms] .= sequences[ms][match_positions[ms]:(match_positions[ms] + motiflen - 1)]\n",
    "        end\n",
    "\n",
    "        # and calculate the consensus score like the paper says - this will become our fitness\n",
    "        fitness_gbest[1], gbest[1] = consensus_score(match_sequences, bg)\n",
    "        \n",
    "        # check if what we got so far is better than in the other runs\n",
    "        if fitness_gbest[1] > final_fitness\n",
    "            final_starts = copy(match_positions)\n",
    "            final_fitness = copy(fitness_gbest[1])\n",
    "            final_motif = StaticArrays.SVector{motiflen, Char}(gbest[1])\n",
    "        end     \n",
    "    end\n",
    "\n",
    "    if includestarts\n",
    "        return(final_fitness, final_motif, final_starts)\n",
    "    else\n",
    "        return(final_fitness, final_motif)\n",
    "    end\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a78195f-2260-4e08-87fe-0129ee8c5b30",
   "metadata": {},
   "source": [
    "# Unit Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c68bcb8-ec4f-4094-955c-7de7598c59bb",
   "metadata": {},
   "source": [
    "## PSO Helper Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80ab44bc-aa95-4bbe-b8bb-3d82413c9c92",
   "metadata": {},
   "source": [
    "### make_profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a241c2d-dcee-4cc0-b536-9ddea5768446",
   "metadata": {},
   "outputs": [],
   "source": [
    "?make_profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7a6e2bd-0974-41fb-852b-6054f305649a",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test make_profile(3, 3, [['A', 'A', 'A'], ['A', 'A', 'A'], ['A', 'A', 'A']] ; pseudo_counts = false) == [1.0 1.0 1.0; 0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0;]\n",
    "    Test.@test make_profile(3, 3, [['C', 'C', 'C'], ['C', 'C', 'C'], ['C', 'C', 'C']] ; pseudo_counts = false) == [0.0 0.0 0.0; 1.0 1.0 1.0; 0.0 0.0 0.0; 0.0 0.0 0.0;]\n",
    "    Test.@test make_profile(3, 3, [['G', 'G', 'G'], ['G', 'G', 'G'], ['G', 'G', 'G']] ; pseudo_counts = false) == [0.0 0.0 0.0; 0.0 0.0 0.0; 1.0 1.0 1.0; 0.0 0.0 0.0;]\n",
    "    Test.@test make_profile(3, 3, [['T', 'T', 'T'], ['T', 'T', 'T'], ['T', 'T', 'T']] ; pseudo_counts = false) == [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0; 1.0 1.0 1.0;]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "163ba920-4588-4d65-9b2f-3d2b6fc13955",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test make_profile(3, 3, [['A', 'A', 'A'], ['A', 'A', 'A'], ['A', 'A', 'A']] ; pseudo_counts = true) == [4/7 4/7 4/7; 1/7 1/7 1/7; 1/7 1/7 1/7; 1/7 1/7 1/7;]\n",
    "    Test.@test make_profile(3, 3, [['C', 'C', 'C'], ['C', 'C', 'C'], ['C', 'C', 'C']] ; pseudo_counts = true) == [1/7 1/7 1/7; 4/7 4/7 4/7; 1/7 1/7 1/7; 1/7 1/7 1/7;]\n",
    "    Test.@test make_profile(3, 3, [['G', 'G', 'G'], ['G', 'G', 'G'], ['G', 'G', 'G']] ; pseudo_counts = true) == [1/7 1/7 1/7; 1/7 1/7 1/7; 4/7 4/7 4/7; 1/7 1/7 1/7;]\n",
    "    Test.@test make_profile(3, 3, [['T', 'T', 'T'], ['T', 'T', 'T'], ['T', 'T', 'T']] ; pseudo_counts = true) == [1/7 1/7 1/7; 1/7 1/7 1/7; 1/7 1/7 1/7; 4/7 4/7 4/7;]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b14ba58-74b4-4fb3-858d-16249aeba806",
   "metadata": {},
   "source": [
    "### background_frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2a026ac-80b3-42e3-b3e0-335ec44e9f3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "?background_frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48e00f1c-6792-4f63-b59f-bc1b6ae3adcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test background_frequency(collect(repeat('A', 100))) == Dict('A' => 100, 'C' =>   0, 'G' =>   0, 'T' =>   0)\n",
    "    Test.@test background_frequency(collect(repeat('C', 100))) == Dict('A' =>   0, 'C' => 100, 'G' =>   0, 'T' =>   0)\n",
    "    Test.@test background_frequency(collect(repeat('G', 100))) == Dict('A' =>   0, 'C' =>   0, 'G' => 100, 'T' =>   0)\n",
    "    Test.@test background_frequency(collect(repeat('T', 100))) == Dict('A' =>   0, 'C' =>   0, 'G' =>   0, 'T' => 100)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e481ce5-5b5b-4ae7-8cf8-d569c3bbb03f",
   "metadata": {},
   "source": [
    "### get_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "286d86d2-cef3-4753-9099-e3c327287018",
   "metadata": {},
   "outputs": [],
   "source": [
    "?get_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9797056b-fee3-468b-b9b6-a03055acf681",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test get_index('A') == 1\n",
    "    Test.@test get_index('C') == 2\n",
    "    Test.@test get_index('G') == 3\n",
    "    Test.@test get_index('T') == 4\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab17e652-56b1-4e3d-90c2-ca617f9402c7",
   "metadata": {},
   "source": [
    "### consensus_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c079a87a-b8ff-4f40-b9cf-615d0b5e7c02",
   "metadata": {},
   "outputs": [],
   "source": [
    "?consensus_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69099332-26fd-4731-89f0-3ba4ebd1a7b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test consensus_score([['A', 'A'], ['A', 'T']], Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)) == (3.0, ['A', 'A'])\n",
    "    Test.@test consensus_score([['A', 'T'], ['A', 'T']], Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)) == (4.0, ['A', 'T'])\n",
    "    Test.@test consensus_score([['T', 'A'], ['A', 'T']], Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)) == (2.0, ['A', 'A'])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5db52f5b-1499-4c92-a3db-a72f3df23dfc",
   "metadata": {
    "tags": []
   },
   "source": [
    "### weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8026c773-93d2-4545-97a2-7d1e36d7e4c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "?weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ceb891a3-152a-45e8-ba70-d0f38fa44dcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test weights(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)) == ones(Float64, 4)\n",
    "    Test.@test weights(Dict('A' => 1/2, 'C' => 1/6, 'G' => 1/6, 'T' => 1/6)) == [1/3, 1, 1, 1]\n",
    "    Test.@test weights(Dict('A' => 1/2, 'C' => 1/2, 'G' => 0, 'T' => 0)) == ones(Float64, 4)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0de3375c-2941-49cc-be81-4dda8adc918c",
   "metadata": {},
   "source": [
    "### update_motif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dad8ba86-9dbe-4d78-9d02-bb448645b0e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "?update_motif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dabd7f7-1463-44be-83d1-e65132d1ef07",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test update_motif(['A', 'A', 'A', 'T'], ['C', 'C', 'C', 'C'], ['G', 'G', 'G', 'G'], [2.0, 1.0, 1.0, 1.0], Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4); deterministic=true) == ['A', 'A', 'A', 'T']\n",
    "    Test.@test update_motif(['A', 'A', 'A', 'T'], ['C', 'C', 'C', 'C'], ['G', 'G', 'G', 'G'], [1.0, 2.0, 1.0, 1.0], Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4); deterministic=true) == ['C', 'C', 'C', 'C']\n",
    "    Test.@test update_motif(['A', 'A', 'A', 'T'], ['C', 'C', 'C', 'C'], ['G', 'G', 'G', 'G'], [1.0, 1.0, 2.0, 1.0], Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4); deterministic=true) == ['G', 'G', 'G', 'G']\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f138841c-0674-47d8-bdc8-d116cf726835",
   "metadata": {},
   "source": [
    "## Our Creative Contributions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f91f2be-4b2d-4e76-9c52-c62ccaf44e8f",
   "metadata": {},
   "source": [
    "### precompute_matchscores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "102ef13a-19c7-4707-8e06-d6094e253cc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "?precompute_matchscores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdbb2cc4-7602-45ee-a9de-7e0ee5a17e48",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)) == Dict(('A', 'A') => 1,\n",
    "        ('A', 'C') => 0,\n",
    "        ('A', 'G') => 0,\n",
    "        ('A', 'T') => 0,\n",
    "        ('C', 'A') => 0,\n",
    "        ('C', 'C') => 1,\n",
    "        ('C', 'G') => 0,\n",
    "        ('C', 'T') => 0,\n",
    "        ('G', 'A') => 0,\n",
    "        ('G', 'C') => 0,\n",
    "        ('G', 'G') => 1,\n",
    "        ('G', 'T') => 0,\n",
    "        ('T', 'A') => 0,\n",
    "        ('T', 'C') => 0,\n",
    "        ('T', 'G') => 0,\n",
    "        ('T', 'T') => 1)\n",
    "        \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1126d96-95db-497c-827d-270eadb6cfed",
   "metadata": {},
   "source": [
    "### precompute_hash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec89ed37-5036-4efc-99f5-a2acff205317",
   "metadata": {},
   "outputs": [],
   "source": [
    "?precompute_hash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9023604f-d5a5-4835-a710-2326d95cc40f",
   "metadata": {},
   "outputs": [],
   "source": [
    "precompute_hash(collect(\"AACGT\"), 3) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6935b598-ff98-47ba-8361-2155f0dae28a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Need to fix this test casee\n",
    "\n",
    "begin\n",
    "    Test.@test precompute_hash(collect(\"AACGT\"), 3)  ==  Dict(2 => Dict('A'=>[1, 0, 0], 'G'=>[0, 0, 1], 'T'=>[0, 0, 0], 'C'=>[0, 1, 0]),\n",
    "        3 => Dict('A'=>[0, 0, 0], 'G'=>[0, 1, 0], 'T'=>[0, 0, 1], 'C'=>[1, 0, 0]),\n",
    "        1 => Dict('A'=>[1, 1, 0], 'G'=>[0, 0, 0], 'T'=>[0, 0, 0], 'C'=>[0, 0, 1]))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35996800-4d95-4170-ad1b-22662142988d",
   "metadata": {},
   "source": [
    "### precompute_motifscores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc5c9f82-8b00-41ed-8f14-efbc29aefc69",
   "metadata": {},
   "outputs": [],
   "source": [
    "?precompute_motifscores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "932693ab-36de-4b91-a7b1-302276e73e5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test precompute_motifscores(['A', 'A', 'A', 'T', 'T', 'T'], precompute_hash(['A', 'A', 'A', 'T', 'T', 'T'], 2), 2, \n",
    "        precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))) ==\n",
    "    Dict(1=> Dict('A' => [1.0, 1.0, 1.0, 0.0, 0.0],\n",
    "    'G' => [0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "    'T' => [0.0, 0.0, 0.0, 1.0, 1.0],\n",
    "    'C' => [0.0, 0.0, 0.0, 0.0, 0.0]),\n",
    "2 => Dict('A' => [1.0, 1.0, 0.0, 0.0, 0.0],\n",
    "    'G' => [0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "    'T' => [0.0, 0.0, 1.0, 1.0, 1.0],\n",
    "    'C' => [0.0, 0.0, 0.0, 0.0, 0.0]))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e1ca56c-4f6c-4157-80fc-2e814f1f8177",
   "metadata": {
    "tags": []
   },
   "source": [
    "### best_match_efficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aed0cfba-6929-42de-82f6-ab9d983a9ed2",
   "metadata": {},
   "outputs": [],
   "source": [
    "?best_match_efficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dcf4122-673a-4928-91a8-cc03c41dfde4",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test best_match_efficient(['A', 'A', 'A', 'T', 'T', 'T'], \n",
    "        precompute_motifscores(['A', 'A', 'A', 'T', 'T', 'T'], precompute_hash(['A', 'A', 'A', 'T', 'T', 'T'], 2), 2, \n",
    "                               precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))), \n",
    "        ['A', 'T'], \n",
    "        precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),\n",
    "        zeros(Float64, 5)) == (2.0, 3)\n",
    "    \n",
    "    Test.@test best_match_efficient(['A', 'A', 'A', 'T', 'T', 'T'], \n",
    "        precompute_motifscores(['A', 'A', 'A', 'T', 'T', 'T'], precompute_hash(['A', 'A', 'A', 'T', 'T', 'T'], 3), 3, \n",
    "                               precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))), \n",
    "        ['A', 'T', 'C'], \n",
    "        precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),\n",
    "        zeros(Float64, 4))  == (2.0, 3)\n",
    "\n",
    "    Test.@test best_match_efficient(['A', 'A', 'A', 'T', 'T', 'T'], \n",
    "        precompute_motifscores(['A', 'A', 'A', 'T', 'T', 'T'], precompute_hash(['A', 'A', 'A', 'T', 'T', 'T'], 4), 4, \n",
    "                               precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))), \n",
    "        ['A', 'T', 'C', 'C'], \n",
    "        precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),\n",
    "        zeros(Float64, 3))  == (2.0, 3)\n",
    "    \n",
    "    Test.@test best_match_efficient(['A', 'A', 'A', 'T', 'T', 'T'], \n",
    "        precompute_motifscores(['A', 'A', 'A', 'T', 'T', 'T'], precompute_hash(['A', 'A', 'A', 'T', 'T', 'T'], 4), 4, \n",
    "                               precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))), \n",
    "        ['A', 'A', 'A', 'T'], \n",
    "        precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),\n",
    "        zeros(Float64, 3))  == (4.0, 1)\n",
    "end\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93b757ca-eabb-4239-8947-aa35d034b3af",
   "metadata": {},
   "source": [
    "### best_possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe1be924-498d-4ee7-ad12-4801e311e81e",
   "metadata": {},
   "outputs": [],
   "source": [
    "?best_possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea2f53ba-ced0-4c56-b42a-f701d1bd23be",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test best_possible([1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4], [2 2; 0 0 ; 0 0 ; 0 0 ;], 2, 10, zeros(Float64, 4, 4), zeros(Float64, 2)) == 4.0\n",
    "    Test.@test best_possible([1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4], [2 2 2; 0 0 0; 0 0 0; 0 0 0;], 2, 10, zeros(Float64, 4, 4), zeros(Float64, 3)) == 6.0\n",
    "    Test.@test best_possible([1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4], [1 1 1; 1 1 1; 0 0 0; 0 0 0;], 2, 10, zeros(Float64, 4, 4), zeros(Float64, 3)) ==\n",
    "    best_possible([1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4], [0 0 0; 0 0 0; 1 1 1; 1 1 1;], 2, 10, zeros(Float64, 4, 4), zeros(Float64, 3))\n",
    "    Test.@test best_possible([1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4], [2 2 2; 1 1 1; 0 0 0; 0 0 0;], 3, 10, zeros(Float64, 4, 4), zeros(Float64, 3)) <\n",
    "    best_possible([1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4; 1/4 1/4 1/4 1/4], [2 2 2; 0 0 0;  0 0 0;  0 0 0;], 2, 10, zeros(Float64, 4, 4), zeros(Float64, 3))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "290a2cad-a0d5-4c4b-a6f0-ded83dd59498",
   "metadata": {},
   "source": [
    "### optimistic_search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9da4b4c9-d895-4b6d-b66a-47f29d152b0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "?optimistic_search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14a4eafc-ff9e-4481-a006-7d7c113bcefa",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    Test.@test optimistic_search([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']],\n",
    "            ['A', 'T', 'T'],\n",
    "            Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4),\n",
    "            -Inf,\n",
    "            precompute_motifscores.([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], \n",
    "            precompute_hash.([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], (3,)), \n",
    "            (3, ), (precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),)),\n",
    "            precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))) == [3, 3, 2]\n",
    "    \n",
    "    Test.@test optimistic_search([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']],\n",
    "            ['A', 'T', 'T'],\n",
    "            Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4),\n",
    "            5.99,\n",
    "            precompute_motifscores.([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], \n",
    "            precompute_hash.([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], (3,)), \n",
    "            (3, ), (precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),)),\n",
    "            precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))) == [3, 3, 2]\n",
    "    \n",
    "    Test.@test optimistic_search([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']],\n",
    "            ['A', 'T', 'T'],\n",
    "            Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4),\n",
    "            6,\n",
    "            precompute_motifscores.([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], \n",
    "            precompute_hash.([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], (3,)), \n",
    "            (3, ), (precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),)),\n",
    "            precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))) == [3, 3, 2]\n",
    "    \n",
    "    Test.@test optimistic_search([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']],\n",
    "            ['A', 'T', 'T'],\n",
    "            Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4),\n",
    "            7,\n",
    "            precompute_motifscores.([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], \n",
    "            precompute_hash.([['A', 'A', 'A', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], (3,)), \n",
    "            (3, ), (precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),)),\n",
    "            precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))) == [3, 0, 0]\n",
    "    \n",
    "    Test.@test optimistic_search([['A', 'A', 'C', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']],\n",
    "            ['A', 'T', 'T'],\n",
    "            Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4),\n",
    "            6,\n",
    "            precompute_motifscores.([['A', 'A', 'C', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], \n",
    "            precompute_hash.([['A', 'A', 'C', 'T', 'T', 'T'], ['C', 'A', 'A', 'T', 'T', 'G'], ['T', 'A', 'T', 'T', 'G', 'G']], (3,)), \n",
    "            (3, ), (precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4)),)),\n",
    "            precompute_matchscores(Dict('A' => 1/4, 'C' => 1/4, 'G' => 1/4, 'T' => 1/4))) == [2, 3, 0]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8187c967-a45f-440b-b32a-94b501c43792",
   "metadata": {},
   "source": [
    "# Benchmarking on Dr. Heber's test data from Piazza"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f4a40f1c-e5cc-45a9-9d99-67f40abe7669",
   "metadata": {},
   "source": [
    "dummy_perfect_k, dummy_perfect_seqs = ReadInputs_class(\"./benchmarking/dummy_data_perfect.txt\")\n",
    "pso_motif(dummy_perfect_seqs, dummy_perfect_k, 5, 50, 10; includestarts = true)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "0b2b3453-5b4e-4e75-9b33-a7bc2ae2b807",
   "metadata": {},
   "source": [
    "dummy_noisy_k, dummy_noisy_seqs = ReadInputs_class(\"./benchmarking/dummy_data_noisy.txt\")\n",
    "pso_motif(dummy_noisy_seqs, dummy_noisy_k, 5, 50, 10; includestarts = true)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebdbfc0a-9ddc-4773-899c-69abc5e88bee",
   "metadata": {},
   "source": [
    "# Benchmarking on Real Data (cotton fiber promoters)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "20ac17ce-60f9-42f1-953c-a92866d089f7",
   "metadata": {},
   "source": [
    "cotton_seqs = ReadInputs(\"./benchmarking/cotton_fiber_promoters.fasta\")"
   ]
  },
  {
   "cell_type": "raw",
   "id": "dd7c68ed-f63f-4d3b-9afa-1162259f0f2e",
   "metadata": {},
   "source": [
    "begin\n",
    "    cotton_k = 6\n",
    "    cotton_score, cotton_motif, cotton_starts = pso_motif(cotton_seqs, cotton_k, 10, 50, 100; includestarts = true)\n",
    "    cotton_match_seqs = [join(cotton_seqs[i][cotton_starts[i]:(cotton_starts[i] + cotton_k - 1)]) for i in 1:length(cotton_seqs)]\n",
    "    println(cotton_score)\n",
    "    println(join(cotton_motif))\n",
    "    println()\n",
    "    for i in cotton_match_seqs\n",
    "        println(i)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "bc06f79f-9d62-42af-8692-338ee9c211f3",
   "metadata": {},
   "source": [
    "# save the sequences to a file\n",
    "open(\"./benchmarking/pso_cotton.txt\", \"w\") do file1\n",
    "    # for printing on the seqLogo website\n",
    "    for i in cotton_match_seqs\n",
    "        println(file1, i)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a506e614-d2c8-41ad-af18-c4fc0eec086f",
   "metadata": {},
   "source": [
    "## compare results to MEME"
   ]
  },
  {
   "cell_type": "raw",
   "id": "9f575fad-73ff-48ce-9b23-d1c313f2d73b",
   "metadata": {},
   "source": [
    "lines = readlines(\"./benchmarking/meme_cotton.txt\")\n",
    "DNA = Vector{Vector{Char}}()\n",
    "for line in lines\n",
    "    push!(DNA, collect(line))\n",
    "end\n",
    "seqs_fromMEME = copy(DNA)\n",
    "\n",
    "lines = readlines(\"./benchmarking/pso_cotton.txt\")\n",
    "DNA = Vector{Vector{Char}}()\n",
    "for line in lines\n",
    "    push!(DNA, collect(line))\n",
    "end\n",
    "seqs_fromPSO = copy(DNA)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "fb8cc8e3-57a9-4180-88eb-85fa7471c5aa",
   "metadata": {},
   "source": [
    "cotton_bg = background_frequency(reduce(hcat, cotton_seqs); prop = true)\n",
    "fit_MEME = consensus_score(seqs_fromMEME, cotton_bg)\n",
    "fit_PSO = consensus_score(seqs_fromPSO, cotton_bg)\n",
    "println(\"MEME fitness: \", fit_MEME, \"; PSO fitness: \", fit_PSO)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28aedc6b-55f6-4be8-a850-55374d38fe8f",
   "metadata": {},
   "source": [
    "# Benchmarking on Synthetic Data"
   ]
  },
  {
   "cell_type": "raw",
   "id": "69d651b9-f750-446f-b231-cf657bdece09",
   "metadata": {},
   "source": [
    "\"\"\"\n",
    "code recycled from gtb, HW2\n",
    "\"\"\"\n",
    "function runTrial(l, d, seq_len, seq_num) \n",
    "    # generate the simulated data\n",
    "    motif, motif_starts, motifs_implanted, sequences = GenerateTestData_ld(seq_num, l, seq_len, d)\n",
    "    \n",
    "    test_score, test_motif = pso_motif(sequences, l, 5, 100, 100)\n",
    "\n",
    "    return(Int(all(test_motif .== motif)))\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "11ff1ff2-0998-43a7-bec9-50a38418d1e0",
   "metadata": {},
   "source": [
    "\"\"\"\n",
    "code recycled from gtb, HW2\n",
    "\"\"\"\n",
    "function benchmarker(l, d, seq_len, seq_num, zz)\n",
    "    println(l, \" \", d, \" \", seq_len, \" \", seq_num)\n",
    "    res = permutedims(reduce(hcat, map(i -> [i for i in @timed(runTrial(l, d, seq_len, seq_num))[(:value, :time, :bytes)]], 1:zz)))\n",
    "    return(StatsBase.mean(res; dims = 1))\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "8f02225e-98d4-4ea5-9334-b5a8b841bce6",
   "metadata": {},
   "source": [
    "# set the parameters of what's going to be tested\n",
    "param_testname = vec(reduce(vcat, [repeat([\"seqlen\"], 4),\n",
    "            repeat([\"seqnum\"], 9),\n",
    "            repeat([\"motiflen\"], 16),\n",
    "            repeat([\"mutations\"], 8)]))\n",
    "param_motiflengths =  vec(reduce(vcat, [repeat([14], 4),\n",
    "            repeat([14], 9),\n",
    "            [i for i in 5:20],\n",
    "            repeat([14], 8)]))\n",
    "param_mutations =  vec(reduce(vcat, [repeat([2], 4),\n",
    "            repeat([2], 9),\n",
    "            repeat([0], 16),\n",
    "            [i for i in 0:7]]))\n",
    "param_seqlenghts =  vec(reduce(vcat, [[100, 500, 1000, 5000],\n",
    "            repeat([500], 9),\n",
    "            repeat([500], 16),\n",
    "            repeat([500], 8)]))\n",
    "param_seqnum =  vec(reduce(vcat, [repeat([100], 4),\n",
    "            [20, 40, 60, 80, 100, 200, 300, 400, 500],\n",
    "            repeat([100], 16),\n",
    "            repeat([100], 8)]))\n",
    "param_numtrials = repeat([100], 4 + 9 + 16 + 8)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "3b540cae-94d3-437f-b4e2-49683eb6915d",
   "metadata": {},
   "source": [
    "# iterate over motif length\n",
    "Random.seed!(15161)\n",
    "bm_res = benchmarker.(param_motiflengths, param_mutations, param_seqlenghts, param_seqnum, param_numtrials)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "e4baecb4-6b79-46f7-82c3-9e36bde5cd11",
   "metadata": {},
   "source": [
    "df_res = DataFrame(reduce(vcat, bm_res), [\"Percent Correct\", \"Time (s)\", \"Memory (Mb)\"])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "e921d3ee-1cb6-46e0-a0f6-be74bea4364f",
   "metadata": {},
   "source": [
    "df_res[!,Symbol(\"Memory (Mb)\")] .= df_res[!,Symbol(\"Memory (Mb)\")] ./ 1000000\n",
    "insertcols!(df_res, 1, Symbol(\"Test Name\") => param_testname)\n",
    "insertcols!(df_res, 1, Symbol(\"Motif Length\") => param_motiflengths)\n",
    "insertcols!(df_res, 1, Symbol(\"Number of Mutations\") => param_mutations)\n",
    "insertcols!(df_res, 1, Symbol(\"Length of Sequences\") => param_seqlenghts)\n",
    "insertcols!(df_res, 1, Symbol(\"Number of Sequences\") => param_seqnum)\n",
    "insertcols!(df_res, 1, Symbol(\"Number of Trials\") => param_numtrials)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b38399b1-2dd5-4e4a-8717-d48b7466b6b2",
   "metadata": {},
   "source": [
    "print(df_res)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f5b5dc50-b2c0-496d-bb2e-734a9ac41a98",
   "metadata": {},
   "source": [
    "# save the data\n",
    "CSV.write(\"./benchmarking/test_results.csv\", df_res)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "4f858117-229b-490a-9fcd-dfeb9bff84e7",
   "metadata": {},
   "source": [
    "df_toplot = CSV.read(\"./benchmarking/test_results.csv\", DataFrame)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45faf81d-f420-44ac-b428-fdabbf79bedc",
   "metadata": {},
   "source": [
    "## Accuracy Plots"
   ]
  },
  {
   "cell_type": "raw",
   "id": "fa0925eb-7c98-41b5-a34a-5367aa749d71",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"seqlen\", df_toplot), \n",
    "    x=Symbol(\"Length of Sequences\"), \n",
    "    y=Symbol(\"Percent Correct\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 1, labels = n -> \"$(round(Int, 100n))%\"), \n",
    "    Gadfly.Guide.xlabel(\"length of DNA sequences (nt)\"),\n",
    "    Gadfly.Guide.ylabel(\"correctly identified\\nmotifs (%)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Accuracy by Sequence Length\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_01.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "e652f706-630c-4e81-aada-bf84f8d10d9f",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"seqnum\", df_toplot), \n",
    "    x=Symbol(\"Number of Sequences\"), \n",
    "    y=Symbol(\"Percent Correct\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 1, labels = n -> \"$(round(Int, 100n))%\"), \n",
    "    Gadfly.Guide.xlabel(\"number of DNA sequences\"),\n",
    "    Gadfly.Guide.ylabel(\"correctly identified\\nmotifs (%)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Accuracy by Number of Sequences\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_02.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "13b99dcd-5c44-4e0f-83e1-516500fe243d",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"motiflen\", df_toplot), \n",
    "    x=Symbol(\"Motif Length\"), \n",
    "    y=Symbol(\"Percent Correct\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 1, labels = n -> \"$(round(Int, 100n))%\"), \n",
    "    Gadfly.Guide.xlabel(\"motif length (nt)\"),\n",
    "    Gadfly.Guide.ylabel(\"correctly identified\\nmotifs (%)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Accuracy by Motif Length\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_03.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f2e67f6e-3907-46a2-abb1-9f9bdf3d19bf",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"mutations\", df_toplot), \n",
    "    x=Symbol(\"Number of Mutations\"), \n",
    "    y=Symbol(\"Percent Correct\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 1, labels = n -> \"$(round(Int, 100n))%\"), \n",
    "    Gadfly.Guide.xlabel(\"number of mutations\\n(x out of 14)\"),\n",
    "    Gadfly.Guide.ylabel(\"correctly identified\\nmotifs (%)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Accuracy by Number of Mutations\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_04.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8727d531-674c-448c-9b32-6d7e18c8120f",
   "metadata": {},
   "source": [
    "## run time/memory plots"
   ]
  },
  {
   "cell_type": "raw",
   "id": "560b5243-d42c-4cad-a784-d02f82aae741",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"seqlen\", df_toplot), \n",
    "    x=Symbol(\"Length of Sequences\"), \n",
    "    y=Symbol(\"Memory (Mb)\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 1200), \n",
    "    Gadfly.Guide.xlabel(\"length of DNA sequences (nt)\"),\n",
    "    Gadfly.Guide.ylabel(\"memory used (Mb)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Memory by Sequence Length\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_05_a.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7ae8320f-c079-457b-8f65-f7c896566ac7",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"seqlen\", df_toplot), \n",
    "    x=Symbol(\"Length of Sequences\"), \n",
    "    y=Symbol(\"Time (s)\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 18), \n",
    "    Gadfly.Guide.xlabel(\"length of DNA sequences (nt)\"),\n",
    "    Gadfly.Guide.ylabel(\"time elapsed (s)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Run Time by Sequence Length\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_05_b.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "e8cf4f19-6651-4108-9d0c-5cd74a1c4a93",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"seqnum\", df_toplot), \n",
    "    x=Symbol(\"Number of Sequences\"),\n",
    "    y=Symbol(\"Memory (Mb)\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 1200), \n",
    "    Gadfly.Guide.xlabel(\"number of DNA sequences\"),\n",
    "    Gadfly.Guide.ylabel(\"memory used (Mb)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Memory by Number of Sequences\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_06_a.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "2f28157f-d5c7-46d5-a10f-75935e96170b",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"seqnum\", df_toplot), \n",
    "    x=Symbol(\"Number of Sequences\"),\n",
    "    y=Symbol(\"Time (s)\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 18), \n",
    "    Gadfly.Guide.xlabel(\"number of DNA sequences\"),\n",
    "    Gadfly.Guide.ylabel(\"time elapsed (s)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Run Time by Number of Sequences\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_06_b.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7a0f72cc-5fdc-4459-99e9-f08b0301a4e3",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"motiflen\", df_toplot), \n",
    "    x=Symbol(\"Motif Length\"), \n",
    "    y=Symbol(\"Memory (Mb)\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 1200), \n",
    "    Gadfly.Guide.xlabel(\"motif length (nt)\"),\n",
    "    Gadfly.Guide.ylabel(\"memory used (Mb)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Memory by Motif Length\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_07_a.svg\", 15cm, 10cm), myplot)\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "50104b2b-d4e0-4403-b614-2aafb1a51c72",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"motiflen\", df_toplot), \n",
    "    x=Symbol(\"Motif Length\"), \n",
    "    y=Symbol(\"Time (s)\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(format = :plain),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 18), \n",
    "    Gadfly.Guide.xlabel(\"motif length (nt)\"),\n",
    "    Gadfly.Guide.ylabel(\"time elapsed (s)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Run Time by Motif Length\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_07_b.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22dc7982-62bc-4468-b658-b9cfa8a6dafc",
   "metadata": {},
   "source": [
    "## Benchmarking for Percent of Sequences Lacking Motif"
   ]
  },
  {
   "cell_type": "raw",
   "id": "973f6744-43c3-4c83-a5a6-2815c0e14750",
   "metadata": {},
   "source": [
    "\"\"\"\n",
    "code recycled from gtb, HW2\n",
    "\"\"\"\n",
    "function runTrial_miss(l, d, seq_len, seq_num, percmiss) \n",
    "    # generate the simulated data\n",
    "    motif, motif_starts, motifs_implanted, sequences = GenerateTestData_ld_somemissing(seq_num, l, seq_len, d, percmiss)\n",
    "    \n",
    "    test_score, test_motif = pso_motif(sequences, l, 5, 100, 100)\n",
    "\n",
    "    return(Int(all(test_motif .== motif)))\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d6b0f35b-0476-4a22-9085-ff9704c1f5f1",
   "metadata": {},
   "source": [
    "\"\"\"\n",
    "code recycled from gtb, HW2\n",
    "\"\"\"\n",
    "function benchmarker_miss(l, d, seq_len, seq_num, zz, pmiss)\n",
    "    println(l, \" \", d, \" \", seq_len, \" \", seq_num, \" \", pmiss)\n",
    "    res = permutedims(reduce(hcat, map(i -> [i for i in @timed(runTrial_miss(l, d, seq_len, seq_num, pmiss))[(:value, :time, :bytes)]], 1:zz)))\n",
    "    return(StatsBase.mean(res; dims = 1))\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "66aa1ad9-4efd-42c5-b569-15f9536eb5f0",
   "metadata": {},
   "source": [
    "# iterate over motif length\n",
    "Random.seed!(84768)\n",
    "bm_miss_res = benchmarker_miss.((14,), (0,), (500,), (100,), (100,), [i/100 for i in 0:10:50])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "6ccb5d97-2e76-4e57-a699-f999c5fece53",
   "metadata": {},
   "source": [
    "Random.seed!(54357)\n",
    "bm_miss_res2 = benchmarker_miss.((14,), (0,), (500,), (100,), (100,), [i/100 for i in 60:10:90])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "bbc8b721-9cd6-417c-8466-176d150c4694",
   "metadata": {},
   "source": [
    "df_miss_res = DataFrame(reduce(vcat, vcat(bm_miss_res, bm_miss_res2)), [\"Percent Correct\", \"Time (s)\", \"Memory (Mb)\"])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "e93adf84-7e95-4590-9d3d-0fa20fd890dc",
   "metadata": {},
   "source": [
    "df_miss_res[!,Symbol(\"Memory (Mb)\")] .= df_miss_res[!,Symbol(\"Memory (Mb)\")] ./ 1000000\n",
    "insertcols!(df_miss_res, 1, Symbol(\"Test Name\") => repeat([\"nomotif\"], 10))\n",
    "insertcols!(df_miss_res, 1, Symbol(\"Motif Length\") => repeat([14], 10))\n",
    "insertcols!(df_miss_res, 1, Symbol(\"Number of Mutations\") => repeat([0], 10))\n",
    "insertcols!(df_miss_res, 1, Symbol(\"Length of Sequences\") => repeat([500], 10))\n",
    "insertcols!(df_miss_res, 1, Symbol(\"Number of Sequences\") => repeat([100], 10))\n",
    "insertcols!(df_miss_res, 1, Symbol(\"Number of Trials\") => repeat([100], 10))\n",
    "insertcols!(df_miss_res, 1, Symbol(\"Percent Sequences Missing Motif\") => [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "4bba5bed-0613-4e93-a6ea-255ee26c14d2",
   "metadata": {},
   "source": [
    "# save the data\n",
    "CSV.write(\"./benchmarking/test_results_nomotif.csv\", df_miss_res)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "66934394-307e-4ddc-bc1b-ace5e26dd61d",
   "metadata": {},
   "source": [
    "df_toplot2 = CSV.read(\"./benchmarking/test_results_nomotif.csv\", DataFrame)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "9b2ad894-e738-4a68-a258-ab7af94e7d8f",
   "metadata": {},
   "source": [
    "Gadfly.set_default_plot_size(15cm, 10cm)\n",
    "\n",
    "myplot = Gadfly.plot(DataFrames.filter(Symbol(\"Test Name\") => n -> n == \"nomotif\", df_toplot2), \n",
    "    x=Symbol(\"Percent Sequences Missing Motif\"), \n",
    "    y=Symbol(\"Percent Correct\"),\n",
    "    Gadfly.Geom.point,\n",
    "    Gadfly.Scale.x_continuous(minvalue = 0, maxvalue = 1, labels = n -> \"$(round(Int, 100n))%\"),\n",
    "    Gadfly.Scale.y_continuous(minvalue = 0, maxvalue = 1, labels = n -> \"$(round(Int, 100n))%\"), \n",
    "    Gadfly.Guide.xlabel(\"sequences with no motif implanted (%)\"),\n",
    "    Gadfly.Guide.ylabel(\"correctly identified\\nmotifs (%)\", orientation =:vertical),\n",
    "    Gadfly.Guide.title(\"Accuracy by Presence of Motif\"),\n",
    "    Gadfly.style(major_label_font_size=18pt, minor_label_font_size=18pt, point_size = 6pt))\n",
    "\n",
    "display(myplot)\n",
    "draw(SVG(\"./figures/Figure_09.svg\", 15cm, 10cm), myplot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d68934d-ace1-45c6-88e5-1a00602ba326",
   "metadata": {},
   "source": [
    "## Compare to MEME\n",
    "Meme commands:\n",
    "```for i in $(ls /mnt/c/Users/grant/OneDrive/Documents/GitHub/CSC530_project/benchmarking/meme_test | sort -V); do /usr/bin/time -ao time.log -f \"%e %M\" ~/meme/bin/meme /mnt/c/Users/grant/OneDrive/Documents/GitHub/CSC530_project/benchmarking/meme_test/\"$i\" -dna -oc . -nostatus -time 14400 -mod oops -nmotifs 1 -minw 14 -maxw 14 -objfun classic -markov_order 0 -text | grep \"^MOTIF \" | awk '{print $2}' >> motifs.txt; done```"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a1a4ab3e-d419-4527-9a05-97b44169fc09",
   "metadata": {},
   "source": [
    "# generat the inputs for meme\n",
    "seq_list = []\n",
    "motif_list = []\n",
    "Random.seed!(423134)\n",
    "for i in 1:100\n",
    "    this_correct_motif, this_motif_starts, this_motifs_implanted, this_sequences = GenerateTestData_ld(100, 14, 500, 2)\n",
    "    append!(seq_list, [this_sequences])\n",
    "    append!(motif_list, [this_correct_motif])\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "32b0c8ac-c283-4378-a424-04ae3b87ed06",
   "metadata": {},
   "source": [
    "# write each of the sequences to a file\n",
    "for i in 1:100\n",
    "    fname = \"./benchmarking/meme_test/test\" * string(i) * \".fasta\"\n",
    "    open(fname, \"w\") do file1\n",
    "        # for printing on the seqLogo website\n",
    "        for j in 1:100\n",
    "            println(file1, \">seq\" * string(j))\n",
    "            println(file1, join(seq_list[i][j]))\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a0781597-ae70-4cd2-b566-524758fa782f",
   "metadata": {},
   "source": [
    "# run the test using our algorithm\n",
    "memevspso_pso = []\n",
    "for i in 1:100\n",
    "    this_outpt1, this_outpt2 = pso_motif(seq_list[i], 14, 5, 100, 100)\n",
    "    append!(memevspso_pso, [this_outpt2])\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "092fdc68-a5c4-40e7-a55a-430229b4ff2e",
   "metadata": {},
   "source": [
    "# accuracu with PSO\n",
    "sum(memevspso_pso .== motif_list) / 100"
   ]
  },
  {
   "cell_type": "raw",
   "id": "9344854c-a148-44cc-a0ad-0104ddee1f0d",
   "metadata": {},
   "source": [
    "# results from PSO\n",
    "# accuracy (%)  time (s)   memory (Mb)\n",
    "#         76%   2.15834      187.272"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d0711267-fbf6-4925-950e-5a797f37dec8",
   "metadata": {},
   "source": [
    "# read in the inputs from MEME\n",
    "memevspso_meme = []\n",
    "for line in readlines(\"C:/Users/grant/OneDrive/Documents/GitHub/CSC530_project/benchmarking/meme_test/motifs.txt\")\n",
    "    append!(memevspso_meme, [collect(line)])\n",
    "end\n",
    "# get the run times\n",
    "memevspso_time = []\n",
    "memevspso_memory = []\n",
    "for line in readlines(\"C:/Users/grant/OneDrive/Documents/GitHub/CSC530_project/benchmarking/meme_test/time.log\")\n",
    "    # read in the parameters from the first line\n",
    "    lines1 = parse.(Float64, split(line))\n",
    "    append!(memevspso_time, lines1[1])\n",
    "    append!(memevspso_memory, lines1[2]/1000)\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f804aed8-f908-451d-82a2-a91c5111b5d2",
   "metadata": {},
   "source": [
    "sum(memevspso_time)/100"
   ]
  },
  {
   "cell_type": "raw",
   "id": "cd34f8d0-1c33-4d15-8a65-7388598fce4d",
   "metadata": {},
   "source": [
    "sum(memevspso_memory)/100"
   ]
  },
  {
   "cell_type": "raw",
   "id": "bc1a2439-39f9-4299-9b21-7ae23118b950",
   "metadata": {},
   "source": [
    "sum(memevspso_meme .== motif_list)/100"
   ]
  },
  {
   "cell_type": "raw",
   "id": "5cfcbf2d-6bf8-4669-b787-ea6492d04bac",
   "metadata": {},
   "source": [
    "# results from PSO\n",
    "# accuracy (%)  time (s)   memory (Mb)\n",
    "#         76%   2.15834      187.272"
   ]
  },
  {
   "cell_type": "raw",
   "id": "556f3e83-f8e3-404d-b5e1-34f240f3cd47",
   "metadata": {},
   "source": [
    "# results from MEME\n",
    "# accuracy (%)  time (s)   memory (Mb)\n",
    "#        100%       4.91         12.67"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6087bf5e-d9cf-41c4-bcb9-d568fa77e8fe",
   "metadata": {},
   "source": [
    "# NOTES TO INSTRUCTORS\n",
    "To run our program, you just need to run the commands below, and change `filetotest = XXXXX` to whatever file you're testing.\n",
    "\n",
    "You need to run all the code chunks above \"5. Unit Testing\" to get these commands to work. `Run All` shoudl also work since the chunks we used for our project report have been silenced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c25b7733-b637-4b6b-8ab4-3c1825c219ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "begin\n",
    "    filetotest = \"./benchmarking/dummy_data_perfect.txt\"\n",
    "    testk, testseqs = ReadInputs_class(filetotest)\n",
    "    x = @timed testscore, testmotif, teststarts = pso_motif(testseqs, testk, 5, 100, 100; includestarts = true)\n",
    "    println(\"Run Time (s)         \", lpad(round(x[:time]; digits= 3), 6, \"0\"))\n",
    "    println(\"Memory Usage (MB)    \", lpad(round(x[:bytes]/1000000; digits= 3), 6, \"0\"))\n",
    "    println()\n",
    "    println(\"Information Content  \", string(testscore))\n",
    "    println(\"Consensus Motif     \", join(testmotif))\n",
    "    println()\n",
    "    println(\"==== Location of Motifs and Their Sequences ====\")\n",
    "    println(\"SeqNum  StartPos    Motif\")\n",
    "    for i in 1:length(teststarts)\n",
    "        println(\"   \", lpad(i, 3, \"0\"), \"      \", lpad(teststarts[i], 4, \"0\"), \"    \", join(testseqs[i][teststarts[i]:(teststarts[i] + length(testmotif) - 1)]))\n",
    "    end\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.2",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.2"
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}