Skip to content

Instantly share code, notes, and snippets.

@abalijepalli
Created February 26, 2015 18:58
Show Gist options
  • Save abalijepalli/f29db46f3c360f2106bf to your computer and use it in GitHub Desktop.
Save abalijepalli/f29db46f3c360f2106bf to your computer and use it in GitHub Desktop.
Plot the results of a MOSAIC (https://github.com/usnistgov/mosaic) analysis.
"""
Plot the results of a MOSAIC analysis.
:Created: 02/26/2015
:Author: Arvind Balijepalli <[email protected]>
:ChangeLog:
.. line-block::
02/26/15 AB Initial version
"""
import mosaic.qdfTrajIO as qdf
import mosaic.abfTrajIO as abf
import mosaic.SingleChannelAnalysis
import mosaic.eventSegment as es
import mosaic.stepResponseAnalysis as sra
import glob
import pylab as pl
import numpy as np
import mosaic.sqlite3MDIO as sql
# Process all ABF files in a directory
mosaic.SingleChannelAnalysis.SingleChannelAnalysis(
'~/ReferenceData/abfSet1',
abf.abfTrajIO,
None,
es.eventSegment,
sra.stepResponseAnalysis
).Run()
# Load the results of the analysis
s=sql.sqlite3MDIO()
s.openDB(glob.glob("~/ReferenceData/abfSet1/*sqlite")[-1])
# We first set up a string that holds the query to retrieve the analysis results. Note that {col}
# will be replaced with the name of the database column when we run the query below.
q = "select {col} from metadata where ProcessingStatus='normal' and ResTime > 0.2 \
and BlockDepth between 0.15 and 0.55"
# Now we run two separate queries - the first returns the blockade depth
# and the second returns the residence time. Note that we simply take the query
# string 'q' above and replace {col} with the column name.
x=np.hstack( s.queryDB( q.format(col='BlockDepth') ) )
y=np.hstack( s.queryDB( q.format(col='ResTime') ) )
# Use matplotlib to plot the results with 2 views:
# i) a 1D histogram of blockade depths and
# ii) a 2D histogram of the residence times vs. blockade depth
fig = pl.gcf()
fig.canvas.set_window_title('Residence Time vs. Blockade Depth')
pl.subplot(2, 1, 1)
pl.hist(x, bins=500, histtype='step', rwidth=0.1)
pl.xticks(())
pl.ylabel("Counts", fontsize=14)
pl.subplot(2, 1, 2)
pl.hist2d(x,y, bins=500)
pl.xlabel("Blockade Depth", fontsize=14)
pl.ylabel("Residence Time (ms)", fontsize=14)
pl.ylim([0.2, 20])
pl.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment