test_PSO_implementation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import Pkg\n",
    "# Pkg.add(\"NBInclude\")\n",
    "# Pkg.add(\"StatsBase\")\n",
    "# Pkg.add(\"Manopt\")\n",
    "# Pkg.add(\"Manifolds\")\n",
    "# Pkg.add(\"Distances\")\n",
    "# import the simulation code\n",
    "import NBInclude\n",
    "import Random\n",
    "import StatsBase\n",
    "#import Manopt\n",
    "import Manifolds\n",
    "import Distances\n",
    "using Manopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the code from the l,d motif simulation notebook\n",
    "NBInclude.@nbinclude(\"simulate_ld_motif.ipynb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the seed\n",
    "Random.seed!(100)\n",
    "# use tuple unpacking to get some test values\n",
    "# NumberOfSequences, LengthMotif, LengthSequences, Distance\n",
    "motif, motif_starts, motifs_implanted, sequences = GenerateTestData_ld(3, 14, 100, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write the scoring function\n",
    "# it is written in such a way that it is continuously defined\n",
    "function Score(M, Starts_inp, Seqs=sequences, Length=length(motif))\n",
    "    maxScore = Length * length(Seqs)\n",
    "    \n",
    "    # handle bad inputs\n",
    "    if Starts_inp isa Vector{Float64} \n",
    "        Starts_down = Int.(floor.(Starts_inp))\n",
    "        Starts_up = Int.(ceil.(Starts_inp))\n",
    "    else\n",
    "        Starts_down = [i for i in Starts_inp]\n",
    "        Starts_up = [i + 1 for i in Starts_inp]\n",
    "    end\n",
    "\n",
    "    Starts_up[Starts_up .< 1] .= 1\n",
    "    Starts_up[Starts_up .> (length(Seqs[1]) - Length + 1)] .= length(Seqs[1]) - Length + 1\n",
    "    \n",
    "    Starts_down[Starts_down .< 1] .= 1\n",
    "    Starts_down[Starts_down .> (length(Seqs[1]) - Length + 1)] .= length(Seqs[1]) - Length + 1\n",
    "\n",
    "    # extract the sequences to iterate over\n",
    "    seqsMatrix_down = permutedims(reduce(hcat, map((s, i) -> s[i:(i + Length - 1)], Seqs, Starts_down)))\n",
    "        \n",
    "    seqsMatrix_up = permutedims(reduce(hcat, map((s, i) -> s[i:(i + Length - 1)], Seqs, Starts_up)))\n",
    "\n",
    "    \n",
    "\n",
    "    \n",
    "    # find the most common element in each column\n",
    "    mostCommon_down = mapslices(StatsBase.mode, seqsMatrix_down, dims = 1)\n",
    "    mostCommon_up = mapslices(StatsBase.mode, seqsMatrix_up, dims = 1)\n",
    "    \n",
    "    # now count how many sequences are equal to the consensus\n",
    "    thisScore_down = sum(map((i, j) -> sum(i .== j), eachslice(seqsMatrix_down, dims = 2), mostCommon_down))\n",
    "    thisScore_up = sum(map((i, j) -> sum(i .== j), eachslice(seqsMatrix_up, dims = 2), mostCommon_up))\n",
    "\n",
    "    \n",
    "    dist1 = Distances.euclidean([i for i in Starts_inp], Starts_down)\n",
    "    dist2 = Distances.euclidean([i for i in Starts_inp], Starts_up)\n",
    "    \n",
    "    thisScore = (1 - (dist1)/(dist1 + dist2)) * thisScore_down + (1 - (dist2)/(dist1 + dist2)) * thisScore_up        \n",
    "\n",
    "    if (isnan(thisScore)) | (thisScore > maxScore)\n",
    "        thisScore = 0\n",
    "    end\n",
    "    # figure out the euclidean distance for Starts from Up and Down\n",
    "    \n",
    "    return(maxScore - thisScore)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iterTup = [1:(len - length(motif) + 1) for len in length.(sequences)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# then use tuple unpacking\n",
    "allScores = map(i -> Score(0, i), Iterators.product(iterTup...))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size(allScores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "findmin(allScores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ok so everything works now\n",
    "# lets try and use particle swarm!\n",
    "M = Manifolds.Euclidean(length(sequences))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# random points\n",
    "n_ = 100\n",
    "pts = map(i -> (length(sequences[1]) - length(motif) - 1) * rand(3) .+ 1, 1:n_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vecs = map(i -> (5 * (2 * rand(3) .- 1)), 1:n_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res =particle_swarm!(M, Score, n = n_, x0 = pts, velocity = vecs, inertia = 0.8, return_options = true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res(M,Score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(Manopt.particle_swarm(M, Score, n = n_, x0 = pts, velocity = vecs, inertia = 0.8, return_options = true))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check how good the answer is\n",
    "Score(0, res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_solver_result(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}