<gh_stars>1-10
#!/usr/bin/python
#
#   HRLAnalysis(TM) Software License - Version 1.0 - August 27th, 2013
#
#   Permission is hereby granted, free of charge, to any person or 
#   organization obtaining a copy of the software and accompanying 
#   documentation covered by this license (the "Software") to use, 
#   reproduce, display, distribute, execute, and transmit the 
#   Software, and to prepare derivative works of the Software, and 
#   to permit third-parties to whom the Software is furnished to do 
#   so, all subject to the following:
#
#   The copyright notices in the Software and this entire statement, 
#   including the above license grant, this restriction and the 
#   following disclaimer, must be included in all copies of the 
#   Software, in whole or in part, and all derivative works of the 
#   Software, unless such copies or derivative works are solely in 
#   the form of machine-executable object code generated by a source 
#   language processor.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
#   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES 
#   OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND 
#   NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR 
#   ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR 
#   OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING 
#   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 
#   OTHER DEALINGS IN THE SOFTWARE, INCLUDING BUT NOT LIMITED TO THE 
#   COMPATIBILITY OF THIS LICENSE WITH OTHER SOFTWARE LICENSES.
#
#
#   runRatAnalysis.py is a script that generates one of the figures 
#   in the HRLAnalysis(TM) paper.  It uses spike data from a rat.  
#   Look at the default data directory (../../test/data/RAT/) for 
#   more details.  This analysis can be ran with no arguments, and produces 
#   a figure ../results/ratFig.eps

import sys
import os
from copy import copy
from optparse import Option, OptionValueError
from optparse import OptionParser

sys.path.append('../lib/')
import libHrlAnalysis

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.lines import Line2D
import matplotlib.transforms as mtransforms
import matplotlib.mpl as mpl

"""
def trace(frame, event, arg):
    print "%s, %s:%d" % (event, frame.f_code.co_filename, frame.f_lineno)
    return trace
"""

def analyzePosAndBurst(options):
    ################ build Analysis ################
    analysis = RatAnalysis(options.startTime, options.endTime, options.spkFile, options.posFile)
    analysis.buildDataStructures()

    ################ plot position ################

    # build axes
    posPlt = plt.axes([.03,.6,.93,.36],aspect='equal')

    # plot 5 ranges of data as partitioned by the timeIntervals.
    # Data within the timeIntervals are examined in more detail, 
    # so denoted by special character / color highlighting 
    tList = options.timeInterval1 + options.timeInterval2
    beg = 0
    idx = beg
    markerStyle = ['o', 'h', 'o', 'h', 'o']
    color = ['w', '#9999CC', 'w', '#99CC99', 'w']
    markerSize = [3.0, 10.0, 3.0, 10.0, 3.0]

    for step in range(4):
        while idx < len(analysis.t_pos) and analysis.t_pos[idx] < tList[step]:
            idx += 1 # find the end of each interval
        # plot the x,y position of the rat
        posPlt.plot(analysis.x_pos[beg:idx+1], analysis.y_pos[beg:idx+1], linestyle='None', color=color[step], marker=markerStyle[step], markersize=markerSize[step], markevery=2, markeredgecolor='none') 
        beg = idx - 1 # overlap the intervals by one data point
    # plot the last interval
    posPlt.plot(analysis.x_pos[beg:], analysis.y_pos[beg:], linestyle='None', color=color[4], marker=markerStyle[4], markersize=markerSize[4], markevery=2, markeredgecolor='none')

    # rat platform was 100cm x 100cm 
    posPlt.set_xlim((-52,52))
    posPlt.set_ylim((-52,52))

    # remove frame
    posPlt.axes.spines['top'].set_visible(False)
    posPlt.axes.spines['left'].set_visible(False)
    posPlt.axes.spines['right'].set_visible(False)
    posPlt.axes.spines['bottom'].set_visible(False)
    posPlt.set_yticks([])
    posPlt.set_xticks([])

    # insert measure bars 
    minX,maxX = posPlt.get_xlim() # should be -52,52
    minY,maxY = posPlt.get_ylim() # should be -52,52

    # x measBar
    l = Line2D([minX-5,minX+15],[minY-5,minY-5],lw=3,color='k') 
    posPlt.add_line(l)
    posPlt.text(minX-2,minY-11, "20 cm",fontsize = 10)
    # y measBar
    l = Line2D([minX-5,minX-5],[minY-5,minY+15],lw=3,color='k') 
    posPlt.add_line(l)
    posPlt.text(minX-23,minY+4, "20 cm",fontsize = 10)

    # reset limits to include measure bars
    posPlt.set_xlim((minX-6,maxX))
    posPlt.set_ylim((minY-6,maxY))



    ################ plot spike raster with bursting ################

    # build axes
    rasterPlt = plt.axes([.03,.35,.93,.22])

    # highlight intervals on which extra analysis is done
    trans = mtransforms.blended_transform_factory(rasterPlt.axes.transData, rasterPlt.axes.transAxes)
    x = np.linspace(options.startTime/1000.0,(options.endTime+1)/1000.0,options.endTime+1-options.startTime)
    fill_intervals = np.zeros((2,len(x)))
    # timeInterval1
    for i in range(options.timeInterval1[0],options.timeInterval1[1]):
        fill_intervals[0][i-options.startTime] = 1
    # the 0 & 20 are lower and upper y-bnds for filling
    rasterPlt.axes.fill_between(x, .1, 20, where=fill_intervals[0]==1, facecolor='#9999CC', transform=trans, edgecolor='none')
    # timeInterval2
    for i in range(options.timeInterval2[0],options.timeInterval2[1]):
        fill_intervals[1][i-options.startTime] = 1
    # the 0 & 20 are lower and upper y-bnds for filling
    rasterPlt.axes.fill_between(x, .1, 20, where=fill_intervals[1]==1, facecolor='#99CC99', transform=trans, edgecolor='none')


    # plot raster
    rasterInfo = analysis.getSpikeTimes()
    times = [t/1000.0 for t in rasterInfo.time]
    rasterPlt.scatter(times, rasterInfo.spikes, color='k',marker='.',linewidth=2)

    # plot bursting
    bursts = analysis.getBursting(options.significance)
    for i, cellBurstInfo in enumerate(bursts):
        for burst in cellBurstInfo:
            burstRange = [burst.start/1000.0, burst.stop/1000.0]
            cellIndex = [i for t in burstRange]
            rasterPlt.plot(burstRange, cellIndex, color='r', alpha=0.7, linewidth=1.5)

    # format
    rasterPlt.axes.set_ylim((analysis.paramsIn().startIdx-.5,analysis.paramsIn().endIdx+.5))
    rasterPlt.axes.set_xlim((options.startTime/1000.0,options.endTime/1000.0))
    rasterPlt.axes.spines['top'].set_visible(False)
    rasterPlt.axes.spines['left'].set_visible(False)
    rasterPlt.axes.spines['right'].set_visible(False)
    rasterPlt.axes.spines['bottom'].set_visible(False)
    rasterPlt.set_yticks([])
    rasterPlt.set_xticks([])

    # x measBar
    minX,maxX = rasterPlt.get_xlim()
    minY,maxY = rasterPlt.get_ylim()
    l = Line2D([minX,minX+2],[minY-.5,minY-.5],lw=3,color='k') 
    rasterPlt.add_line(l)
    rasterPlt.text(minX+.3,minY-2, "2 s",fontsize = 10)
    rasterPlt.set_ylim(minY-2,maxY)

def analyzePairSynch((startTime, endTime), options, horzOffset, border_color):
    # build colored border to match plot with analyzed intervals
    border_pairSynchPlt = plt.axes([horzOffset,.03,.43,.29],aspect='equal')
    border_pairSynchPlt.axes.set_ylim(0,10)
    border_pairSynchPlt.axes.set_xlim(0,10)
    x = np.linspace(0,10,1000)
    trans = mtransforms.blended_transform_factory(border_pairSynchPlt.axes.transData, border_pairSynchPlt.axes.transAxes)
    border_pairSynchPlt.axes.spines['top'].set_visible(False)
    border_pairSynchPlt.axes.spines['left'].set_visible(False)
    border_pairSynchPlt.axes.spines['right'].set_visible(False)
    border_pairSynchPlt.axes.spines['bottom'].set_visible(False)
    border_pairSynchPlt.set_yticks([])
    border_pairSynchPlt.set_xticks([])

    # the 0 & 20 are lower and upper bounds for filling
    border_pairSynchPlt.axes.fill_between(x, 0, 20, facecolor=border_color, transform=trans, edgecolor='none')

    # build analysis structures
    analysis = RatAnalysis(startTime, endTime, options.spkFile, options.posFile)
    analysis.buildDataStructures()
    analysis.setEmptyTrainSynchVal(1.0)
    
    # construct axes
    pairSynchPlt = plt.axes([horzOffset+.01,.03+.01,.43-.02,.29-.02],aspect='equal')

    # plot pairwise neuron synchrony for each neuron pair for current time interval
    numNeurons = analysis.paramsIn().endIdx - analysis.paramsIn().startIdx
    synchronyGrid = np.zeros((numNeurons,numNeurons))

    # Get the synchrony value for each neuron pair (i,j).
    # 0 is considered synchronous, where 1 is the opposite
    for i in range(numNeurons):
        for j in range(numNeurons):
            if i == j:
                synchronyGrid[i][j] = 0.0
            else:
                synchronyInfo = analysis.getPairSynchrony(analysis.paramsIn().startIdx + i, analysis.paramsIn().startIdx + j)
                val = analysis.calcSPIKEDistanceAvg(synchronyInfo)
                synchronyGrid[i][j] = val

    # normalize the color scheme to give us the colors we want
    norm = mpl.colors.Normalize(vmin=-.3, vmax=1.0)
    pairSynchPlt.pcolormesh(synchronyGrid, cmap=plt.get_cmap('afmhot'), norm=norm)    

    # remove frame
    pairSynchPlt.axes.spines['top'].set_visible(False)
    pairSynchPlt.axes.spines['left'].set_visible(False)
    pairSynchPlt.axes.spines['right'].set_visible(False)
    pairSynchPlt.axes.spines['bottom'].set_visible(False)
    pairSynchPlt.set_yticks([])
    pairSynchPlt.set_xticks([])



class RatAnalysis(libHrlAnalysis.HrlNeuralAnalysis):
    def __init__(self, startTime, endTime, spkFile, posFile):
        libHrlAnalysis.HrlNeuralAnalysis.__init__(self,0,0,0,0,libHrlAnalysis.vector_string())
        self.paramsIn().startTime = startTime
        self.paramsIn().endTime = endTime
        self.spkFile = spkFile
        self.posFile = posFile
        self.x_pos = []
        self.y_pos = []
        self.t_pos = []

    def buildDataStructures(self):
        # add spiketrains
        for cell,spTrain in enumerate(open(self.spkFile)):
            self.cellActivity().append(libHrlAnalysis.vector_int())
            # put data in [ms]
            for spikeT in [int(round(eval(i)*1000)) for i in spTrain.split()]:
                if self.paramsIn().startTime <= spikeT and spikeT <= self.paramsIn().endTime:
                    self.cellActivity()[cell].append(spikeT)
                    self.spikeActivity().append(libHrlAnalysis.int_pair(spikeT,cell))
        self.sortSpikeActivity()
        
        # generate position data
        for t,x,y in [eval(i) for i in open(self.posFile)]:
            time = int(round(t*1000))
            if self.paramsIn().startTime <= time and time <= self.paramsIn().endTime:
                self.x_pos.append(x)
                self.y_pos.append(y)
                self.t_pos.append(time)

        self.paramsIn().startIdx = 0
        self.paramsIn().endIdx = cell
        self.paramsIn().isDataCompiled = True
        return True

def check_tuple(option, opt, value):
    try:
        tup = eval(value)
        if len(tup) == 2 and tup[0] > 0 and tup[0] < tup[1] and tup[1] < 600000:
            return tuple(tup)
        else:
            raise OptionValueError("option %s: invalid time-range tuple: %r" % (opt, value))
    except:
        raise OptionValueError("option %s: invalid time-range tuple: %r" % (opt, value))

class MyOption (Option):
    TYPES = Option.TYPES + ("tuple",)
    TYPE_CHECKER = copy(Option.TYPE_CHECKER)
    TYPE_CHECKER["tuple"] = check_tuple

def getParser():
    parser = OptionParser(option_class=MyOption)
            
    parser.add_option("-o", "--outputPath", dest="outputPath",
                      action="store",type="string",default="../results/",
                      help="Output path to store the plots.", 
                      metavar="PATH")
    
    parser.add_option("-s", "--startTime", dest="startTime",
                      action="store",type="int",default=411500,
                      help="The start time of the analysis [ms], range: [0..595000].", 
                      metavar="TYPE")

    parser.add_option("-e", "--endTime", dest="endTime",
                      action="store",type="int",default=466000,
                      help="The end time of the analysis [ms], range: [5..600000].", 
                      metavar="TYPE")

    parser.add_option("-d", "--spikeFile", dest="spkFile",
                      action="store",type="string",default='../../test/data/RAT/data-06020601.txt',
                      help="The file containing the spike data.", 
                      metavar="FILE")

    parser.add_option("-p", "--posFile", dest="posFile",
                      action="store",type="string",default='../../test/data/RAT/pos-06020601.txt',
                      help="The file containing the rat position data.", 
                      metavar="FILE")

    parser.add_option("-A", "--analysisIntervalA", dest="timeInterval1",
                      action="store",type="tuple",default=(422500,425000),
                      help="First analysis interval [ms]: (begTime, endTime).", 
                      metavar="TYPE")

    parser.add_option("-B", "--analysisIntervalB", dest="timeInterval2",
                      action="store",type="tuple",default=(462500,465000),
                      help="Second analysis interval [ms]: (begTime, endTime).", 
                      metavar="TYPE")
    
    parser.add_option("-b", "--burstSig", dest="significance",
                      action="store",type="float",default=.05,
                      help="Bursting significance.",
                      metavar="TYPE")
    
    

    return parser

if __name__ == "__main__":
    #sys.settrace(trace)

    # Parse the command line.
    parser = getParser()
    (options,args) = parser.parse_args()

    # build figure frame
    fig = plt.figure(1, figsize=(2*3.34646,2*4.13386))

    # add labels
    fig.text(0.25,0.96, "A.",ha='center',fontsize = 14)
    fig.text(0.02,0.575, "B.",ha='center',fontsize = 14)
    fig.text(0.02,0.32, "C.",ha='center',fontsize = 14)

    
    horzOffset1 = .05
    horzOffset2 = .52
    
    analyzePosAndBurst(options)
    analyzePairSynch(options.timeInterval1, options, horzOffset1, '#9999CC')
    analyzePairSynch(options.timeInterval2, options, horzOffset2, '#99CC99')

    # save figure
    plt.savefig(options.outputPath + 'ratFig.eps', dpi=1200)


#!/usr/bin/python
#
#   HRLAnalysis(TM) Software License - Version 1.0 - August 27th, 2013
#
#   Permission is hereby granted, free of charge, to any person or 
#   organization obtaining a copy of the software and accompanying 
#   documentation covered by this license (the "Software") to use, 
#   reproduce, display, distribute, execute, and transmit the 
#   Software, and to prepare derivative works of the Software, and 
#   to permit third-parties to whom the Software is furnished to do 
#   so, all subject to the following:
#
#   The copyright notices in the Software and this entire statement, 
#   including the above license grant, this restriction and the 
#   following disclaimer, must be included in all copies of the 
#   Software, in whole or in part, and all derivative works of the 
#   Software, unless such copies or derivative works are solely in 
#   the form of machine-executable object code generated by a source 
#   language processor.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
#   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES 
#   OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND 
#   NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR 
#   ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR 
#   OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING 
#   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 
#   OTHER DEALINGS IN THE SOFTWARE, INCLUDING BUT NOT LIMITED TO THE 
#   COMPATIBILITY OF THIS LICENSE WITH OTHER SOFTWARE LICENSES.
#
#
#   runRatAnalysis.py is a script that generates one of the figures 
#   in the HRLAnalysis(TM) paper.  It uses spike data from a rat.  
#   Look at the default data directory (../../test/data/RAT/) for 
#   more details.  This analysis can be ran with no arguments, and produces 
#   a figure ../results/ratFig.eps

import sys
import os
from copy import copy
from optparse import Option, OptionValueError
from optparse import OptionParser

sys.path.append('../lib/')
import libHrlAnalysis

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.lines import Line2D
import matplotlib.transforms as mtransforms
import matplotlib.mpl as mpl

"""
def trace(frame, event, arg):
    print "%s, %s:%d" % (event, frame.f_code.co_filename, frame.f_lineno)
    return trace
"""

def analyzePosAndBurst(options):
    ################ build Analysis ################
    analysis = RatAnalysis(options.startTime, options.endTime, options.spkFile, options.posFile)
    analysis.buildDataStructures()

    ################ plot position ################

    # build axes
    posPlt = plt.axes([.03,.6,.93,.36],aspect='equal')

    # plot 5 ranges of data as partitioned by the timeIntervals.
    # Data within the timeIntervals are examined in more detail, 
    # so denoted by special character / color highlighting 
    tList = options.timeInterval1 + options.timeInterval2
    beg = 0
    idx = beg
    markerStyle = ['o', 'h', 'o', 'h', 'o']
    color = ['w', '#9999CC', 'w', '#99CC99', 'w']
    markerSize = [3.0, 10.0, 3.0, 10.0, 3.0]

    for step in range(4):
        while idx < len(analysis.t_pos) and analysis.t_pos[idx] < tList[step]:
            idx += 1 # find the end of each interval
        # plot the x,y position of the rat
        posPlt.plot(analysis.x_pos[beg:idx+1], analysis.y_pos[beg:idx+1], linestyle='None', color=color[step], marker=markerStyle[step], markersize=markerSize[step], markevery=2, markeredgecolor='none') 
        beg = idx - 1 # overlap the intervals by one data point
    # plot the last interval
    posPlt.plot(analysis.x_pos[beg:], analysis.y_pos[beg:], linestyle='None', color=color[4], marker=markerStyle[4], markersize=markerSize[4], markevery=2, markeredgecolor='none')

    # rat platform was 100cm x 100cm 
    posPlt.set_xlim((-52,52))
    posPlt.set_ylim((-52,52))

    # remove frame
    posPlt.axes.spines['top'].set_visible(False)
    posPlt.axes.spines['left'].set_visible(False)
    posPlt.axes.spines['right'].set_visible(False)
    posPlt.axes.spines['bottom'].set_visible(False)
    posPlt.set_yticks([])
    posPlt.set_xticks([])

    # insert measure bars 
    minX,maxX = posPlt.get_xlim() # should be -52,52
    minY,maxY = posPlt.get_ylim() # should be -52,52

    # x measBar
    l = Line2D([minX-5,minX+15],[minY-5,minY-5],lw=3,color='k') 
    posPlt.add_line(l)
    posPlt.text(minX-2,minY-11, "20 cm",fontsize = 10)
    # y measBar
    l = Line2D([minX-5,minX-5],[minY-5,minY+15],lw=3,color='k') 
    posPlt.add_line(l)
    posPlt.text(minX-23,minY+4, "20 cm",fontsize = 10)

    # reset limits to include measure bars
    posPlt.set_xlim((minX-6,maxX))
    posPlt.set_ylim((minY-6,maxY))



    ################ plot spike raster with bursting ################

    # build axes
    rasterPlt = plt.axes([.03,.35,.93,.22])

    # highlight intervals on which extra analysis is done
    trans = mtransforms.blended_transform_factory(rasterPlt.axes.transData, rasterPlt.axes.transAxes)
    x = np.linspace(options.startTime/1000.0,(options.endTime+1)/1000.0,options.endTime+1-options.startTime)
    fill_intervals = np.zeros((2,len(x)))
    # timeInterval1
    for i in range(options.timeInterval1[0],options.timeInterval1[1]):
        fill_intervals[0][i-options.startTime] = 1
    # the 0 & 20 are lower and upper y-bnds for filling
    rasterPlt.axes.fill_between(x, .1, 20, where=fill_intervals[0]==1, facecolor='#9999CC', transform=trans, edgecolor='none')
    # timeInterval2
    for i in range(options.timeInterval2[0],options.timeInterval2[1]):
        fill_intervals[1][i-options.startTime] = 1
    # the 0 & 20 are lower and upper y-bnds for filling
    rasterPlt.axes.fill_between(x, .1, 20, where=fill_intervals[1]==1, facecolor='#99CC99', transform=trans, edgecolor='none')


    # plot raster
    rasterInfo = analysis.getSpikeTimes()
    times = [t/1000.0 for t in rasterInfo.time]
    rasterPlt.scatter(times, rasterInfo.spikes, color='k',marker='.',linewidth=2)

    # plot bursting
    bursts = analysis.getBursting(options.significance)
    for i, cellBurstInfo in enumerate(bursts):
        for burst in cellBurstInfo:
            burstRange = [burst.start/1000.0, burst.stop/1000.0]
            cellIndex = [i for t in burstRange]
            rasterPlt.plot(burstRange, cellIndex, color='r', alpha=0.7, linewidth=1.5)

    # format
    rasterPlt.axes.set_ylim((analysis.paramsIn().startIdx-.5,analysis.paramsIn().endIdx+.5))
    rasterPlt.axes.set_xlim((options.startTime/1000.0,options.endTime/1000.0))
    rasterPlt.axes.spines['top'].set_visible(False)
    rasterPlt.axes.spines['left'].set_visible(False)
    rasterPlt.axes.spines['right'].set_visible(False)
    rasterPlt.axes.spines['bottom'].set_visible(False)
    rasterPlt.set_yticks([])
    rasterPlt.set_xticks([])

    # x measBar
    minX,maxX = rasterPlt.get_xlim()
    minY,maxY = rasterPlt.get_ylim()
    l = Line2D([minX,minX+2],[minY-.5,minY-.5],lw=3,color='k') 
    rasterPlt.add_line(l)
    rasterPlt.text(minX+.3,minY-2, "2 s",fontsize = 10)
    rasterPlt.set_ylim(minY-2,maxY)

def analyzePairSynch((startTime, endTime), options, horzOffset, border_color):
    # build colored border to match plot with analyzed intervals
    border_pairSynchPlt = plt.axes([horzOffset,.03,.43,.29],aspect='equal')
    border_pairSynchPlt.axes.set_ylim(0,10)
    border_pairSynchPlt.axes.set_xlim(0,10)
    x = np.linspace(0,10,1000)
    trans = mtransforms.blended_transform_factory(border_pairSynchPlt.axes.transData, border_pairSynchPlt.axes.transAxes)
    border_pairSynchPlt.axes.spines['top'].set_visible(False)
    border_pairSynchPlt.axes.spines['left'].set_visible(False)
    border_pairSynchPlt.axes.spines['right'].set_visible(False)
    border_pairSynchPlt.axes.spines['bottom'].set_visible(False)
    border_pairSynchPlt.set_yticks([])
    border_pairSynchPlt.set_xticks([])

    # the 0 & 20 are lower and upper bounds for filling
    border_pairSynchPlt.axes.fill_between(x, 0, 20, facecolor=border_color, transform=trans, edgecolor='none')

    # build analysis structures
    analysis = RatAnalysis(startTime, endTime, options.spkFile, options.posFile)
    analysis.buildDataStructures()
    analysis.setEmptyTrainSynchVal(1.0)
    
    # construct axes
    pairSynchPlt = plt.axes([horzOffset+.01,.03+.01,.43-.02,.29-.02],aspect='equal')

    # plot pairwise neuron synchrony for each neuron pair for current time interval
    numNeurons = analysis.paramsIn().endIdx - analysis.paramsIn().startIdx
    synchronyGrid = np.zeros((numNeurons,numNeurons))

    # Get the synchrony value for each neuron pair (i,j).
    # 0 is considered synchronous, where 1 is the opposite
    for i in range(numNeurons):
        for j in range(numNeurons):
            if i == j:
                synchronyGrid[i][j] = 0.0
            else:
                synchronyInfo = analysis.getPairSynchrony(analysis.paramsIn().startIdx + i, analysis.paramsIn().startIdx + j)
                val = analysis.calcSPIKEDistanceAvg(synchronyInfo)
                synchronyGrid[i][j] = val

    # normalize the color scheme to give us the colors we want
    norm = mpl.colors.Normalize(vmin=-.3, vmax=1.0)
    pairSynchPlt.pcolormesh(synchronyGrid, cmap=plt.get_cmap('afmhot'), norm=norm)    

    # remove frame
    pairSynchPlt.axes.spines['top'].set_visible(False)
    pairSynchPlt.axes.spines['left'].set_visible(False)
    pairSynchPlt.axes.spines['right'].set_visible(False)
    pairSynchPlt.axes.spines['bottom'].set_visible(False)
    pairSynchPlt.set_yticks([])
    pairSynchPlt.set_xticks([])



class RatAnalysis(libHrlAnalysis.HrlNeuralAnalysis):
    def __init__(self, startTime, endTime, spkFile, posFile):
        libHrlAnalysis.HrlNeuralAnalysis.__init__(self,0,0,0,0,libHrlAnalysis.vector_string())
        self.paramsIn().startTime = startTime
        self.paramsIn().endTime = endTime
        self.spkFile = spkFile
        self.posFile = posFile
        self.x_pos = []
        self.y_pos = []
        self.t_pos = []

    def buildDataStructures(self):
        # add spiketrains
        for cell,spTrain in enumerate(open(self.spkFile)):
            self.cellActivity().append(libHrlAnalysis.vector_int())
            # put data in [ms]
            for spikeT in [int(round(eval(i)*1000)) for i in spTrain.split()]:
                if self.paramsIn().startTime <= spikeT and spikeT <= self.paramsIn().endTime:
                    self.cellActivity()[cell].append(spikeT)
                    self.spikeActivity().append(libHrlAnalysis.int_pair(spikeT,cell))
        self.sortSpikeActivity()
        
        # generate position data
        for t,x,y in [eval(i) for i in open(self.posFile)]:
            time = int(round(t*1000))
            if self.paramsIn().startTime <= time and time <= self.paramsIn().endTime:
                self.x_pos.append(x)
                self.y_pos.append(y)
                self.t_pos.append(time)

        self.paramsIn().startIdx = 0
        self.paramsIn().endIdx = cell
        self.paramsIn().isDataCompiled = True
        return True

def check_tuple(option, opt, value):
    try:
        tup = eval(value)
        if len(tup) == 2 and tup[0] > 0 and tup[0] < tup[1] and tup[1] < 600000:
            return tuple(tup)
        else:
            raise OptionValueError("option %s: invalid time-range tuple: %r" % (opt, value))
    except:
        raise OptionValueError("option %s: invalid time-range tuple: %r" % (opt, value))

class MyOption (Option):
    TYPES = Option.TYPES + ("tuple",)
    TYPE_CHECKER = copy(Option.TYPE_CHECKER)
    TYPE_CHECKER["tuple"] = check_tuple

def getParser():
    parser = OptionParser(option_class=MyOption)
            
    parser.add_option("-o", "--outputPath", dest="outputPath",
                      action="store",type="string",default="../results/",
                      help="Output path to store the plots.", 
                      metavar="PATH")
    
    parser.add_option("-s", "--startTime", dest="startTime",
                      action="store",type="int",default=411500,
                      help="The start time of the analysis [ms], range: [0..595000].", 
                      metavar="TYPE")

    parser.add_option("-e", "--endTime", dest="endTime",
                      action="store",type="int",default=466000,
                      help="The end time of the analysis [ms], range: [5..600000].", 
                      metavar="TYPE")

    parser.add_option("-d", "--spikeFile", dest="spkFile",
                      action="store",type="string",default='../../test/data/RAT/data-06020601.txt',
                      help="The file containing the spike data.", 
                      metavar="FILE")

    parser.add_option("-p", "--posFile", dest="posFile",
                      action="store",type="string",default='../../test/data/RAT/pos-06020601.txt',
                      help="The file containing the rat position data.", 
                      metavar="FILE")

    parser.add_option("-A", "--analysisIntervalA", dest="timeInterval1",
                      action="store",type="tuple",default=(422500,425000),
                      help="First analysis interval [ms]: (begTime, endTime).", 
                      metavar="TYPE")

    parser.add_option("-B", "--analysisIntervalB", dest="timeInterval2",
                      action="store",type="tuple",default=(462500,465000),
                      help="Second analysis interval [ms]: (begTime, endTime).", 
                      metavar="TYPE")
    
    parser.add_option("-b", "--burstSig", dest="significance",
                      action="store",type="float",default=.05,
                      help="Bursting significance.",
                      metavar="TYPE")
    
    

    return parser

if __name__ == "__main__":
    #sys.settrace(trace)

    # Parse the command line.
    parser = getParser()
    (options,args) = parser.parse_args()

    # build figure frame
    fig = plt.figure(1, figsize=(2*3.34646,2*4.13386))

    # add labels
    fig.text(0.25,0.96, "A.",ha='center',fontsize = 14)
    fig.text(0.02,0.575, "B.",ha='center',fontsize = 14)
    fig.text(0.02,0.32, "C.",ha='center',fontsize = 14)

    
    horzOffset1 = .05
    horzOffset2 = .52
    
    analyzePosAndBurst(options)
    analyzePairSynch(options.timeInterval1, options, horzOffset1, '#9999CC')
    analyzePairSynch(options.timeInterval2, options, horzOffset2, '#99CC99')

    # save figure
    plt.savefig(options.outputPath + 'ratFig.eps', dpi=1200)
#!/usr/bin/python # #   HRLAnalysis(TM) Software License - Version 1.0 - August 27th, 2013 # #   Permission is hereby granted, free of charge, to any person or #   organization obtaining a copy of the software and accompanying #   documentation covered by this license (the "Software") to use, #   reproduce, display, distribute, execute, and transmit the #   Software, and to prepare derivative works of the Software, and #   to permit third-parties to whom the Software is furnished to do #   so, all subject to the following: # #   The copyright notices in the Software and this entire statement, #   including the above license grant, this restriction and the #   following disclaimer, must be included in all copies of the #   Software, in whole or in part, and all derivative works of the #   Software, unless such copies or derivative works are solely in #   the form of machine-executable object code generated by a source #   language processor. # #   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, #   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES #   OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND #   NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR #   ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR #   OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING #   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR #   OTHER DEALINGS IN THE SOFTWARE, INCLUDING BUT NOT LIMITED TO THE #   COMPATIBILITY OF THIS LICENSE WITH OTHER SOFTWARE LICENSES. # # #   runRatAnalysis.py is a script that generates one of the figures #   in the HRLAnalysis(TM) paper.  It uses spike data from a rat. #   Look at the default data directory (../../test/data/RAT/) for #   more details.  This analysis can be ran with no arguments, and produces #   a figure ../results/ratFig.eps def trace(frame, event, arg):     print "%s, %s:%d" % (event, frame.f_code.co_filename, frame.f_lineno)     return trace ################ build Analysis ################ ################ plot position ################ # build axes # plot 5 ranges of data as partitioned by the timeIntervals. # Data within the timeIntervals are examined in more detail, # so denoted by special character / color highlighting # find the end of each interval # plot the x,y position of the rat # overlap the intervals by one data point # plot the last interval # rat platform was 100cm x 100cm # remove frame # insert measure bars # should be -52,52 # should be -52,52 # x measBar # y measBar # reset limits to include measure bars ################ plot spike raster with bursting ################ # build axes # highlight intervals on which extra analysis is done # timeInterval1 # the 0 & 20 are lower and upper y-bnds for filling # timeInterval2 # the 0 & 20 are lower and upper y-bnds for filling # plot raster # plot bursting # format # x measBar # build colored border to match plot with analyzed intervals # the 0 & 20 are lower and upper bounds for filling # build analysis structures # construct axes # plot pairwise neuron synchrony for each neuron pair for current time interval # Get the synchrony value for each neuron pair (i,j). # 0 is considered synchronous, where 1 is the opposite # normalize the color scheme to give us the colors we want # remove frame # add spiketrains # put data in [ms] # generate position data #sys.settrace(trace) # Parse the command line. # build figure frame # add labels # save figure