Monday, May 14, 2018

Dynamical systems model of cardio-respiratory interactions

Currently I'm working on a model of cardio-respiratory interactions. Given human experimental data, the model reproduces the average heart rate, the average respiratory phase duration (inspiration/expiration), and the Respiratory Sinus Arrhythmia (RSA). The most sophisticated part of the model is the respiratory central pattern generator (previously published). My current goal is to optimize the model so that it reproduces RSA for several individuals.

After I get it working with the current dataset-- which was measured while participants watched Disney's Fantasia-- the next step is to fit the model to real clinical data from septic ER patients. The end goal is for our collaborators to use our model as a basis for a machine learning algorithm that predicts the risk and most likely cause of death for septic patients.

Below are some pretty pictures generated by the model.

Monday, January 16, 2017

Haiku Twitter Bot

Over summer 2016, I decided to make a Twitter bot just for fun. Here, I'll describe the bot's programming at a high level. But first, a little background...

I started undergraduate as an english major... then graduated with a degree in neuroscience. Although I have pretty firmly switched my career plans toward math and science, I still appreciate and respect art. As a reflection on my artistic pipedream, I designed my bot to write haikus. If you don't remember haikus from high school, hit that link in the previous sentence.

For connecting to Twitter, I used the tweepy library. I simply wrote a subclass of StreamListener to process incoming tweets and write them to a file called tweets.txt. You can read more about this process in the tweepy docs (linked above).

Each time I start my bot, most of the high-level functionality is located in one function, do_tweet().

For counting syllables and determining parts of speech, I used NLTK.

Counting syllables:


Choosing words for the haiku:

As you can see, the algorithm randomly chooses a word of syllable length sw, which is generated by the pick_syl() function. If this is hard to interpret, don't worry. It took me a while to come up with an algorithm that would randomly choose words but still adhere to the 5/7/5 haiku format. The conditions like if ix == 0 and line == 0 are there to determine which line of the haiku is being written. In this case, the first word of the first line is being chosen. Then the lines:

sdict1 = dict([(sk, sv) for sk, sv in sdict.items() if pos_tag(word_tokenize(sk),
tagset='universal')[0][1]=='NOUN'])
# Chooses a noun
gsyls = syl_of_size(sdict1, sw) # Chooses a noun of syllable length sw

Then, setting up the haiku in order and writing to the file haiku.txt:

So currently this bot would write haikus with randomly chosen words, which is cool... but we can do better! ;^)
Eventually when I have a little more time, I'll give meta ai ai haiku more smarts, but for now the bot uses a simple algorithm for determining the "best" tweets, then recycles words from these tweets. Dumb, I know... it's a work in progress.


The function get_best() chooses tweets with the highest "weight" to include in the bot's corpus for writing future haikus.

Here's a link to my Twitter bot meta ai ai haiku: robcapps.com/docs/haiku. Please feel free to ask questions in the comments section!

Monday, September 12, 2016

Using plt.imshow to label timeseries data with matplotlib

Recently I've been trying to find the most effective way to color-code time series data using matplotlib.

For instance, if I have a raster plot like this:
...and I want to label it like this:
How?

After lots of failed attempts involving horrific for loops and ax.fill_between, ax.fill, and several other inferior solutions, I found a quick and easy fix:

After plotting your time series raster plot,

    Given: 
      - an array colors in (R, G, B) format (same length as your time series data)
      - nr_neurons is the number of axes (in this case, 4)
      - vspace is a float representing the space between each time axis
      - A timeline variable, tspace (created with np.linspace)


Note the parameter aspect='auto'. This makes sure the image takes the shape of your existing axis.

Tuesday, August 23, 2016

QGIS Area of Availability Plugin

8/23/16

I'm working on a set of QGIS plugins, and one that I think should be out there for grabs is an area of availability calculator. (EDIT: you should be able to install the plugin via the public QGIS repository soon)

Given a start point and a maximum distance to travel, we want to find the boundary of accessible regions via a road network.

Boundary of area of availability, plotted as points


We build a road network graph and calculate the road distances using Dijkstra's algorithm. We can then find the points corresponding to the outer boundary of the area of availability.

The code for this calculation is below (taken from this page):

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

pStart = QgsPoint(65.5462, 57.1509)
delta = qgis.utils.iface.mapCanvas().getCoordinateTransform().mapUnitsPerPixel() * 1

rb = QgsRubberBand(qgis.utils.iface.mapCanvas(), True)
rb.setColor(Qt.green)
rb.addPoint(QgsPoint(pStart.x() - delta, pStart.y() - delta))
rb.addPoint(QgsPoint(pStart.x() + delta, pStart.y() - delta))
rb.addPoint(QgsPoint(pStart.x() + delta, pStart.y() + delta))
rb.addPoint(QgsPoint(pStart.x() - delta, pStart.y() + delta))

tiedPoints = director.makeGraph(builder, [pStart])
graph = builder.graph()
tStart = tiedPoints[0]

idStart = graph.findVertex(tStart)

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)

upperBound = []
r = 2000.0
i = 0
while i < len(cost):
  if cost[i] > r and tree[i] != -1:
    outVertexId = graph.arc(tree [i]).outVertex()
    if cost[outVertexId] < r:
      upperBound.append(i)
  i = i + 1

for i in upperBound:
  centerPoint = graph.vertex(i).point()
  rb = QgsRubberBand(qgis.utils.iface.mapCanvas(), True)
  rb.setColor(Qt.red)
  rb.addPoint(QgsPoint(centerPoint.x() - delta, centerPoint.y() - delta))
  rb.addPoint(QgsPoint(centerPoint.x() + delta, centerPoint.y() - delta))
  rb.addPoint(QgsPoint(centerPoint.x() + delta, centerPoint.y() + delta))
  rb.addPoint(QgsPoint(centerPoint.x() - delta, centerPoint.y() + delta))

I'm not quite finished with my plugin implementation, but you can check out the source code on my GitHub.

If you want to try your hand at developing a QGIS plugin, I strongly recommend starting with Plugin Builder.

Sunday, July 31, 2016

7/31/16, Current Work
  Since December 2015, I've been collaborating on a research publication. The publication will be centered on a Python software package called FreqPy.
  FreqPy is a tool designed to assist with qualitative analysis of neural firing rate data and to bolster arguments for neural population coding in subpopulations of neurons.
  We're currently in the verification/justification stages; I synthesized some sine waves as test data to verify that our methodology is working and just finished creating a set of figures.
  Check out the video below: