| 6 | | import numpy |
|---|
| 7 | | import os |
|---|
| 8 | | try : |
|---|
| 9 | | import pylab |
|---|
| 10 | | except Exception: |
|---|
| 11 | | print "Warning: pylab not present" |
|---|
| 12 | | |
|---|
| 13 | | |
|---|
| 14 | | class SpikeTrain(object): |
|---|
| 15 | | """ |
|---|
| 16 | | This class defines a spike train as a list of the events times. |
|---|
| 17 | | |
|---|
| 18 | | Event times are given in a list (sparse representation) in milliseconds. |
|---|
| 19 | | |
|---|
| 20 | | >>> s1 = SpikeTrain([0.0, 0.1, 0.2, 0.5]) |
|---|
| 21 | | >>> s2 = SpikeTrain(numpy.array([0, 1, 2, 5]), dt=0.1) |
|---|
| 22 | | >>> assert all(s1.spike_times == s2.spike_times) |
|---|
| 23 | | >>> print s1.format(relative=True) |
|---|
| 24 | | [ 0. 0.1 0.1 0.3] |
|---|
| 25 | | >>> s1.isi() |
|---|
| 26 | | array([ 0.1, 0.1, 0.3]) |
|---|
| 27 | | >>> s1.mean_rate() |
|---|
| 28 | | 8.0 |
|---|
| 29 | | >>> s1.cv_isi() |
|---|
| 30 | | 0.565685424949 |
|---|
| 31 | | """ |
|---|
| 32 | | # TODO in the definition, should a spike train be ordered? |
|---|
| 33 | | |
|---|
| 34 | | # TODO : should we handle different population shapes differently? |
|---|
| 35 | | |
|---|
| 36 | | ####################################################################### |
|---|
| 37 | | ## Constructor and key methods to manipulate the SpikeTrain objects ## |
|---|
| 38 | | ####################################################################### |
|---|
| 39 | | def __init__(self, spike_times, dt=None, t_start=None, t_stop=None): |
|---|
| 40 | | """ |
|---|
| 41 | | `spike_times` is a list/numpy array of spike times (in milliseconds) |
|---|
| 42 | | |
|---|
| 43 | | TODO: We proposed initially that : |
|---|
| 44 | | If `dt`is specified, the time values should be ints, |
|---|
| 45 | | and will be multiplied by `dt`, otherwise time values should be floats. |
|---|
| 46 | | Markus suggested that the internal representation should be integer and that analysis methods also. |
|---|
| 47 | | |
|---|
| 48 | | If `t_start` and `t_stop` are not specified, they are inferred from the data. |
|---|
| 49 | | |
|---|
| 50 | | the format adopted for the representation is relative=False, quantized=False |
|---|
| 51 | | |
|---|
| 52 | | """ |
|---|
| 53 | | |
|---|
| 54 | | if dt is not None: |
|---|
| 55 | | if dt<=0: |
|---|
| 56 | | raise ValueError("dt must be greater than zero") |
|---|
| 57 | | #self.spike_times *= dt |
|---|
| 58 | | |
|---|
| 59 | | self.dt = dt |
|---|
| 60 | | self.t_start = t_start |
|---|
| 61 | | self.t_stop = t_stop |
|---|
| 62 | | |
|---|
| 63 | | self.spike_times = numpy.array(spike_times, 'float') |
|---|
| 64 | | |
|---|
| 65 | | # If t_start is not None, we resize the spike_train keeping only |
|---|
| 66 | | # the spikes with t >= t_start |
|---|
| 67 | | if self.t_start is not None: |
|---|
| 68 | | idx = numpy.where(self.spike_times >= self.t_start)[0] |
|---|
| 69 | | self.spike_times = self.spike_times[idx] |
|---|
| 70 | | assert numpy.all(self.spike_times >= self.t_start), \ |
|---|
| 71 | | "Spike times out of range (t_start=%s, min(spike_times)=%s" % (self.t_start,self.spike_times.min()) |
|---|
| 72 | | |
|---|
| 73 | | # If t_stop is not None, we resize the spike_train keeping only |
|---|
| 74 | | # the spikes with t <= t_stop |
|---|
| 75 | | if self.t_stop is not None: |
|---|
| 76 | | idx = numpy.where(self.spike_times <= self.t_stop)[0] |
|---|
| 77 | | self.spike_times = self.spike_times[idx] |
|---|
| 78 | | assert numpy.all(self.spike_times <= self.t_stop), \ |
|---|
| 79 | | "Spike times out of range (t_stop=%s, max(spike_times)=%s" % (self.t_stop,self.spike_times.max()) |
|---|
| 80 | | |
|---|
| 81 | | # Here we deal with the t_start and t_stop values if the SpikeTrain |
|---|
| 82 | | # is empty, with only one element or several elements, if we |
|---|
| 83 | | # need to guess t_start and t_stop |
|---|
| 84 | | # no element : t_start = 0, t_stop = dt if dt is defined, else 1.0 ms |
|---|
| 85 | | # 1 element : t_start = time, t_stop = time + dt |
|---|
| 86 | | # several : t_start = min(time), t_stop = max(time) |
|---|
| 87 | | |
|---|
| 88 | | size = len(spike_times) |
|---|
| 89 | | if size == 0: |
|---|
| 90 | | if self.t_start is None: |
|---|
| 91 | | self.t_start = 0 |
|---|
| 92 | | if self.t_stop is None: |
|---|
| 93 | | self.t_stop = self.dt or 1.0 |
|---|
| 94 | | elif size == 1: # spike list may be empty |
|---|
| 95 | | if self.t_start is None: |
|---|
| 96 | | self.t_start = self.spike_times[0] |
|---|
| 97 | | if self.t_stop is None: |
|---|
| 98 | | self.t_stop = self.spike_times[0] + (self.dt or 1.0) |
|---|
| 99 | | elif size > 1: |
|---|
| 100 | | if self.t_start is None: |
|---|
| 101 | | self.t_start = numpy.min(self.spike_times) |
|---|
| 102 | | if numpy.any(self.spike_times < self.t_start): |
|---|
| 103 | | raise ValueError("Spike times must not be less than t_start") |
|---|
| 104 | | if self.t_stop is None: |
|---|
| 105 | | self.t_stop = numpy.max(self.spike_times) |
|---|
| 106 | | if numpy.any(self.spike_times > self.t_stop): |
|---|
| 107 | | raise ValueError("Spike times must not be greater than t_stop") |
|---|
| 108 | | |
|---|
| 109 | | if self.t_start >= self.t_stop : |
|---|
| 110 | | raise Exception("Incompatible time interval for the creation of the SpikeTrain. t_start=%s, t_stop=%s" % (self.t_start, self.t_stop)) |
|---|
| 111 | | if self.t_start < 0: |
|---|
| 112 | | raise ValueError("t_start must not be negative") |
|---|
| 113 | | if numpy.any(self.spike_times < 0): |
|---|
| 114 | | raise ValueError("Spike times must not be negative") |
|---|
| 115 | | |
|---|
| 116 | | def __str__(self): |
|---|
| 117 | | return str(self.spike_times) |
|---|
| 118 | | |
|---|
| 119 | | def __len__(self): |
|---|
| 120 | | return len(self.spike_times) |
|---|
| 121 | | |
|---|
| 122 | | def format(self, relative=False, quantized=False): |
|---|
| 123 | | """ |
|---|
| 124 | | a function to format a spike train from one format to another. |
|---|
| 125 | | outputs a numpy array |
|---|
| 126 | | |
|---|
| 127 | | """ |
|---|
| 128 | | spike_times = self.spike_times.copy() |
|---|
| 129 | | spike_times.sort() |
|---|
| 130 | | |
|---|
| 131 | | if relative and len(spike_times)>0: |
|---|
| 132 | | spike_times[1:] = spike_times[1:] - spike_times[:-1] |
|---|
| 133 | | |
|---|
| 134 | | if quantized: |
|---|
| 135 | | assert quantized > 0, "quantized must either be False or a positive number" |
|---|
| 136 | | #spike_times = numpy.array([time/self.quantized for time in spike_times],int) |
|---|
| 137 | | spike_times = (spike_times/quantized).round().astype('int') |
|---|
| 138 | | |
|---|
| 139 | | return spike_times |
|---|
| 140 | | |
|---|
| 141 | | ####################################################################### |
|---|
| 142 | | ## Analysis methods that can be applied to a SpikeTrain object ## |
|---|
| 143 | | ####################################################################### |
|---|
| 144 | | def isi(self): |
|---|
| 145 | | # TODO this needs some thinking to know how to handle the border, in particular the |
|---|
| 146 | | # first spike and t_start |
|---|
| 147 | | return self.format(relative=True, quantized=False)[1:] |
|---|
| 148 | | |
|---|
| 149 | | # Returns the mean firing rate of the SpikeTrain |
|---|
| 150 | | def mean_rate(self, t_start=None, t_stop=None): |
|---|
| 151 | | """ Mean firing rate between t_start and t_stop in Hz |
|---|
| 152 | | |
|---|
| 153 | | NOTE: avoided calling where when defaults settings t_start=self.t_start, t_stop=self.t_stop |
|---|
| 154 | | """ |
|---|
| 155 | | if (t_start==None) & (t_stop==None): |
|---|
| 156 | | t_start=self.t_start |
|---|
| 157 | | t_stop=self.t_stop |
|---|
| 158 | | idx = self.spike_times |
|---|
| 159 | | else: |
|---|
| 160 | | if t_start==None: t_start=self.t_start |
|---|
| 161 | | if t_stop==None: t_stop=self.t_stop |
|---|
| 162 | | idx = numpy.where((self.spike_times >= t_start) & (self.spike_times <= t_stop))[0] |
|---|
| 163 | | return 1000.*len(idx)/(t_stop-t_start) |
|---|
| 164 | | |
|---|
| 165 | | # Returns the cv of the isi of the SpikeTrain |
|---|
| 166 | | def cv_isi(self): |
|---|
| 167 | | """ |
|---|
| 168 | | cv_isi is the ratio between the standard deviation and the mean of the ISI |
|---|
| 169 | | |
|---|
| 170 | | The irregularity of individual spike trains is measured by the squared |
|---|
| 171 | | coefficient of variation of the corresponding inter-spike interval (ISI) |
|---|
| 172 | | distribution normalized by the square of its mean. |
|---|
| 173 | | In point processes, low values reflect more regular spiking, a |
|---|
| 174 | | clock-like pattern yields CV2= 0. On the other hand, CV2 = 1 indicates |
|---|
| 175 | | Poisson-type behavior. As a measure for irregularity in the network one |
|---|
| 176 | | can use the average irregularity across all neurons. |
|---|
| 177 | | |
|---|
| 178 | | TODO: is it useful to get the std of CV? |
|---|
| 179 | | """ |
|---|
| 180 | | isi = self.isi() |
|---|
| 181 | | #TODO return an exception if mean.isi= 0 |
|---|
| 182 | | if len(isi) > 0: |
|---|
| 183 | | return numpy.std(isi)/numpy.mean(isi) |
|---|
| 184 | | else: |
|---|
| 185 | | raise Exception("No spikes in the SpikeTrain !") |
|---|
| 186 | | |
|---|
| 187 | | def fano_factor_isi(self): |
|---|
| 188 | | """ returns the fano factor of this spike trains ISI (see SpikeList.fano_factor)""" |
|---|
| 189 | | isi = self.isi() |
|---|
| 190 | | if len(isi) > 0: |
|---|
| 191 | | #firing_rate = self.time_histogram(time_bin,False) |
|---|
| 192 | | fano = numpy.var(isi)/numpy.mean(isi) |
|---|
| 193 | | return fano |
|---|
| 194 | | else: |
|---|
| 195 | | raise Exception("No spikes in the SpikeTrain !") |
|---|
| 196 | | |
|---|
| 197 | | def time_axis(self, time_bin): |
|---|
| 198 | | """ |
|---|
| 199 | | Return a time axis between t_start and t_stop according to a time_bin |
|---|
| 200 | | """ |
|---|
| 201 | | return numpy.arange(self.t_start, self.t_stop, time_bin) |
|---|
| 202 | | |
|---|
| 203 | | def raster_plot(self, t_start=None, t_stop=None, color='b'): |
|---|
| 204 | | """ |
|---|
| 205 | | Generate a raster plot with the SpikeTrain in a subwindow of interest, |
|---|
| 206 | | defined by t_start and t_stop. Those are not the global t_start and t_stop |
|---|
| 207 | | of the SpikeTrain objects. If not defined, we use the one of the SpikeTrain |
|---|
| 208 | | object |
|---|
| 209 | | """ |
|---|
| 210 | | if t_start is None: |
|---|
| 211 | | t_start = self.t_start |
|---|
| 212 | | if t_stop is None: |
|---|
| 213 | | t_stop = self.t_stop |
|---|
| 214 | | idx = numpy.where((self.spike_times >= t_start) & (self.spike_times <= t_stop))[0] |
|---|
| 215 | | spikes = self.spike_times[idx] |
|---|
| 216 | | if len(spikes) > 0: |
|---|
| 217 | | pylab.figure() |
|---|
| 218 | | pylab.scatter(spikes,numpy.ones(len(spikes)), c=color) |
|---|
| 219 | | pylab.xlabel("Time (ms)", size="x-large") |
|---|
| 220 | | pylab.ylabel("Neuron #", size="x-large") |
|---|
| 221 | | |
|---|
| 222 | | # Method to sort a SpikeTrain |
|---|
| 223 | | def sort_by_time(self): |
|---|
| 224 | | self.spike_time = sort(self.spike_time) |
|---|
| 225 | | |
|---|
| 226 | | def subSpikeTrain(self, t_start, t_stop): |
|---|
| 227 | | """ Returns a spike train sliced between t_start and t_stop |
|---|
| 228 | | t_start and t_stop may either be single values or sequences of start |
|---|
| 229 | | and stop times. |
|---|
| 230 | | """ |
|---|
| 231 | | if hasattr(t_start, '__len__'): |
|---|
| 232 | | assert len(t_start) == len(t_stop) |
|---|
| 233 | | mask = False |
|---|
| 234 | | for t0,t1 in zip(t_start, t_stop): |
|---|
| 235 | | mask = mask | (self.spike_times >= t0) & (self.spike_times <= t1) |
|---|
| 236 | | t_start = t_start[0] |
|---|
| 237 | | t_stop = t_stop[-1] |
|---|
| 238 | | else: |
|---|
| 239 | | mask = (self.spike_times >= t_start) & (self.spike_times <= t_stop) |
|---|
| 240 | | idx = numpy.where(mask)[0] |
|---|
| 241 | | spikes = self.spike_times[idx] |
|---|
| 242 | | return SpikeTrain(spikes, self.dt, t_start, t_stop) |
|---|
| 243 | | |
|---|
| 244 | | def time_histogram(self, time_bin , normalized=True): |
|---|
| 245 | | """ |
|---|
| 246 | | Bin the spikes with the specified bin width. The first and last bins |
|---|
| 247 | | are calculated from `self.t_start` and `self.t_stop`. |
|---|
| 248 | | If `normalized` is True, the bin values are scaled so as to represent |
|---|
| 249 | | firing rates in spikes/second, otherwise it's the number of spikes per bin. |
|---|
| 250 | | """ |
|---|
| 251 | | nbins = self.time_axis(time_bin) |
|---|
| 252 | | (hist, lower_edges ) = numpy.histogram(self.spike_times, nbins) |
|---|
| 253 | | if normalized: |
|---|
| 254 | | hist *= 1000.0/time_bin |
|---|
| 255 | | return hist |
|---|
| 256 | | |
|---|
| 257 | | def rescale(self): |
|---|
| 258 | | """Rescale spike times to make t_start = 0.""" |
|---|
| 259 | | if self.t_start != 0: |
|---|
| 260 | | self.spike_times -= self.t_start |
|---|
| 261 | | self.t_stop -= self.t_start |
|---|
| 262 | | self.t_start = 0.0 |
|---|
| 263 | | |
|---|
| 264 | | def tuning_curve(self, var_array, normalized=False, method='sum'): |
|---|
| 265 | | """ |
|---|
| 266 | | Calculate a firing-rate tuning curve with respect to some variable. |
|---|
| 267 | | Assumes that some variable, such as stimulus orientation, varies |
|---|
| 268 | | throughout the recording. The values taken by this variable should be |
|---|
| 269 | | supplied in a numpy array `var_array`. The spike train is binned |
|---|
| 270 | | according to the number of values in `var_array`, e.g., if there are |
|---|
| 271 | | N values, the spikes are binned with a bin width |
|---|
| 272 | | (`self.t_stop`-`self.t_start`)/N |
|---|
| 273 | | so that each bin is associated with a particular value of the variable |
|---|
| 274 | | in `var_array`. |
|---|
| 275 | | |
|---|
| 276 | | The return value is a dictionary whose keys are the distinct values in |
|---|
| 277 | | `val_array`. The values in the dict depend on the arguments `method` and |
|---|
| 278 | | `normalized`. |
|---|
| 279 | | |
|---|
| 280 | | If `normalized` is False, the responses (bin values) are spike counts, |
|---|
| 281 | | if True, they are firing rates. |
|---|
| 282 | | If `method` == "max", each value is the maximum response for a given |
|---|
| 283 | | value of the variable. |
|---|
| 284 | | If `method` == "sum", each value is the summed response... |
|---|
| 285 | | If `method` == "mean", ... you get the idea. |
|---|
| 286 | | |
|---|
| 287 | | (If someone can rewrite this more clearly, please do so!) |
|---|
| 288 | | """ |
|---|
| 289 | | binwidth = (self.t_stop - self.t_start)/len(var_array) |
|---|
| 290 | | time_histogram = self.time_histogram(binwidth, normalized=normalized) |
|---|
| 291 | | assert len(time_histogram) == len(var_array) |
|---|
| 292 | | tuning_curve = {} |
|---|
| 293 | | counts = {} |
|---|
| 294 | | for k, x in zip(var_array, time_histogram): |
|---|
| 295 | | if not tuning_curve.has_key(k): |
|---|
| 296 | | tuning_curve[k] = 0 |
|---|
| 297 | | counts[k] = 0 |
|---|
| 298 | | if method in ('sum', 'mean'): |
|---|
| 299 | | tuning_curve[k] += x |
|---|
| 300 | | counts[k] += 1 |
|---|
| 301 | | elif method == 'max': |
|---|
| 302 | | tuning_curve[k] = max(x, tuning_curve[k]) |
|---|
| 303 | | else: |
|---|
| 304 | | raise Exception() |
|---|
| 305 | | if method == 'mean': |
|---|
| 306 | | for k in tuning_curve.keys(): |
|---|
| 307 | | tuning_curve[k] = tuning_curve[k]/counts[k] |
|---|
| 308 | | return tuning_curve |
|---|
| 309 | | |
|---|
| 310 | | |
|---|
| 311 | | class SpikeList(object): |
|---|
| 312 | | """ |
|---|
| 313 | | A SpikeList object is a list of SpikeTrain objects. |
|---|
| 314 | | |
|---|
| 315 | | >>> sl = SpikeList(3, [(0, 0.1), (1, 0.1), (0, 0.2)]) |
|---|
| 316 | | >>> type( sl.spiketrains[0] ) |
|---|
| 317 | | <type SpikeTrain> |
|---|
| 318 | | >>> sl.spiketrains[0].spike_times |
|---|
| 319 | | array([ 0.1, 0.2]) |
|---|
| 320 | | >>> sl.as_ids_times() |
|---|
| 321 | | (array([0,0,1]), array([0.1,0.2,0.1])) |
|---|
| 322 | | >>> sl.as_ids_times(relative=True) |
|---|
| 323 | | (array([0,1,0]), array([0.1,0.0,0.1])) |
|---|
| 324 | | >>> sl.as_ids_times(quantized=0.01) |
|---|
| 325 | | (array([0,0,1]), array([10,20,10])) |
|---|
| 326 | | >>> sl.as_list_id_time() |
|---|
| 327 | | [(0,0.1), (0,0.2), (1,0.1)] |
|---|
| 328 | | >>> sl.as_id_list_times() |
|---|
| 329 | | [(0, array([0.1, 0.2])), (1, array([0.1]))] |
|---|
| 330 | | >>> sl.as_time_list_ids() |
|---|
| 331 | | [(0.1, [0,1]), (0.2, [0])] |
|---|
| 332 | | >>> sl.as_2byN_array() |
|---|
| 333 | | >>> array([[0,0,1],[0.1,0.2,0.1]]) |
|---|
| 334 | | |
|---|
| 335 | | """ |
|---|
| 336 | | ####################################################################### |
|---|
| 337 | | ## Constructor and key methods to manipulate the SpikeList objects ## |
|---|
| 338 | | ####################################################################### |
|---|
| 339 | | def __init__(self, spikes, id_list, dt=None, t_start=None, t_stop=None, label=''): |
|---|
| 340 | | """ |
|---|
| 341 | | `spikes` is a list/tuple of (id,t) tuples (id in id_list) |
|---|
| 342 | | `id_list` is the list of ids of all recorded cells (needed for silent cells) |
|---|
| 343 | | If `dt`is specified, the time values should be ints, |
|---|
| 344 | | and will be multiplied by `dt`, otherwise time values should be floats. |
|---|
| 345 | | If `t_start` and `t_stop` are not specified, they are inferred from the data. |
|---|
| 346 | | |
|---|
| 347 | | dt, t_start and t_stop are shared for all SpikeTrains |
|---|
| 348 | | |
|---|
| 349 | | """ |
|---|
| 350 | | self.id_list = id_list |
|---|
| 351 | | self.t_start = t_start |
|---|
| 352 | | self.t_stop = t_stop |
|---|
| 353 | | self.dt = dt |
|---|
| 354 | | self.label = label |
|---|
| 355 | | # transform spikes in a spike array |
|---|
| 356 | | self.spiketrains = {} |
|---|
| 357 | | for id in id_list: |
|---|
| 358 | | self.spiketrains[id] = [] |
|---|
| 359 | | for id,time in spikes: |
|---|
| 360 | | if id in id_list: #id_list can be a subset of the list of recorded neurons |
|---|
| 361 | | self.spiketrains[id].append(time) |
|---|
| 362 | | |
|---|
| 363 | | # writing as a list of SpikeTrains |
|---|
| 364 | | for id,spikes in self.spiketrains.items(): # |
|---|
| 365 | | self.spiketrains[id] = SpikeTrain(spike_times=spikes, dt=self.dt, t_start=self.t_start, t_stop=self.t_stop) |
|---|
| 366 | | if len(self) > 0 and (self.t_start is None or self.t_stop is None): |
|---|
| 367 | | self.__calc_startstop() |
|---|
| 368 | | |
|---|
| 369 | | def N(self): |
|---|
| 370 | | return len(self.id_list) |
|---|
| 371 | | |
|---|
| 372 | | def __calc_startstop(self): |
|---|
| 373 | | """ |
|---|
| 374 | | t_start and t_stop are shared for all neurons, so we take min and max values respectively. |
|---|
| 375 | | TO DO : check the t_start and t_stop parameters for a SpikeList. Is it commun to |
|---|
| 376 | | all the spikeTrains within the spikelist or each spikelistes do need its own. |
|---|
| 377 | | """ |
|---|
| 378 | | if len(self) > 0: |
|---|
| 379 | | if self.t_start is None: |
|---|
| 380 | | start_times = numpy.array([self.spiketrains[idx].t_start for idx in self.id_list]) |
|---|
| 381 | | self.t_start = start_times.min() |
|---|
| 382 | | for id in self.spiketrains.keys(): |
|---|
| 383 | | self.spiketrains[id].t_start = self.t_start |
|---|
| 384 | | if self.t_stop is None: |
|---|
| 385 | | stop_times = numpy.array([self.spiketrains[idx].t_stop for idx in self.id_list]) |
|---|
| 386 | | self.t_stop = stop_times.max() |
|---|
| 387 | | for id in self.spiketrains.keys(): |
|---|
| 388 | | self.spiketrains[id].t_stop = self.t_stop |
|---|
| 389 | | else: |
|---|
| 390 | | raise Exception("No SpikeTrains") |
|---|
| 391 | | |
|---|
| 392 | | def __getitem__(self, i): |
|---|
| 393 | | return self.spiketrains[i] |
|---|
| 394 | | |
|---|
| 395 | | def __setitem__(self, i, val): |
|---|
| 396 | | assert isinstance(val, SpikeTrain), "A SpikeList object can only contain SpikeTrain objects" |
|---|
| 397 | | self.spiketrains[i] = val |
|---|
| 398 | | self.id_list.append(i) |
|---|
| 399 | | self.__calc_startstop() |
|---|
| 400 | | |
|---|
| 401 | | def __len__(self): |
|---|
| 402 | | return len(self.spiketrains) |
|---|
| 403 | | |
|---|
| 404 | | def append(self, id, spiketrain): |
|---|
| 405 | | """ |
|---|
| 406 | | Add a SpikeTrain object to the SpikeList |
|---|
| 407 | | """ |
|---|
| 408 | | assert isinstance(spiketrain, SpikeTrain), "A SpikeList object can only contain SpikeTrain objects" |
|---|
| 409 | | if id in self.id_list: |
|---|
| 410 | | raise Exception("Id already present in SpikeList.Use setitem instead()") |
|---|
| 411 | | else: |
|---|
| 412 | | self.spiketrains[id] = spiketrain |
|---|
| 413 | | self.id_list.append(id) |
|---|
| 414 | | self.N += 1 |
|---|
| 415 | | self.__calc_startstop() |
|---|
| 416 | | |
|---|
| 417 | | def get_time_parameters(self): |
|---|
| 418 | | """ |
|---|
| 419 | | Returns the time parameters of the SpikeList (t_start, t_stop, dt) |
|---|
| 420 | | """ |
|---|
| 421 | | return (self.t_start, self.t_stop, self.dt) |
|---|
| 422 | | |
|---|
| 423 | | # Same as for the SpikeTrain object |
|---|
| 424 | | def time_axis(self, time_bin): |
|---|
| 425 | | return numpy.arange(self.t_start, self.t_stop, time_bin) |
|---|
| 426 | | |
|---|
| 427 | | def concatenate(self, SpikeList_list): |
|---|
| 428 | | """ |
|---|
| 429 | | Concatenation of a list of SpikeLists to the current SpikeList |
|---|
| 430 | | """ |
|---|
| 431 | | # We check that Spike Lists have similar time_axis |
|---|
| 432 | | sl_= SpikeList_list[0] |
|---|
| 433 | | for sl in SpikeList_list: |
|---|
| 434 | | if not sl.get_time_parameters() == self.get_time_parameters(): |
|---|
| 435 | | raise Exception("Spike Lists should have similar time_axis") |
|---|
| 436 | | for sl in SpikeList_list: |
|---|
| 437 | | for id in sl.id_list: |
|---|
| 438 | | self.append(id, sl.spiketrains[id]) |
|---|
| 439 | | |
|---|
| 440 | | |
|---|
| 441 | | def idsubSpikeList(self, id_list): |
|---|
| 442 | | """ |
|---|
| 443 | | Generate a new SpikeList truncated from a particular sublist of cells |
|---|
| 444 | | """ |
|---|
| 445 | | # We check what are the elements that are in self.id_list and not in |
|---|
| 446 | | # id_list. We remove such elements from the SpikeList |
|---|
| 447 | | new_SpkList = SpikeList([], [], self.dt, self.t_start, self.t_stop) |
|---|
| 448 | | if isinstance(id_list,int): |
|---|
| 449 | | id_list = numpy.random.permutation(self.id_list)[0:id_list] |
|---|
| 450 | | for id in id_list: |
|---|
| 451 | | try: |
|---|
| 452 | | new_SpkList.append(id, self.spiketrains[id]) |
|---|
| 453 | | except Exception: |
|---|
| 454 | | print "Item %s is not in the source SpikeList" %id |
|---|
| 455 | | return new_SpkList |
|---|
| 456 | | |
|---|
| 457 | | def timesubSpikeList(self, t_start, t_stop): |
|---|
| 458 | | """ |
|---|
| 459 | | Generate a new SpikeList truncated from particular time boundaries |
|---|
| 460 | | returns a new SpikeList |
|---|
| 461 | | |
|---|
| 462 | | """ |
|---|
| 463 | | new_SpkList = SpikeList([], [], self.dt, t_start, t_stop) |
|---|
| 464 | | for id in self.id_list: |
|---|
| 465 | | new_SpkList.append(id, self.spiketrains[id].subSpikeTrain(t_start, t_stop)) |
|---|
| 466 | | new_SpkList.__calc_startstop() |
|---|
| 467 | | return new_SpkList |
|---|
| 468 | | |
|---|
| 469 | | |
|---|
| 470 | | ####################################################################### |
|---|
| 471 | | ## Analysis methods that can be applied to a SpikeTrain object ## |
|---|
| 472 | | ####################################################################### |
|---|
| 473 | | |
|---|
| 474 | | |
|---|
| 475 | | def isi(self, nbins=100, display=False): |
|---|
| 476 | | """ |
|---|
| 477 | | Return the list of all the isi vectors for all the SpikeTrains objects |
|---|
| 478 | | within the SpikeList. If display is True, then it plots the distribution |
|---|
| 479 | | of the ISI |
|---|
| 480 | | """ |
|---|
| 481 | | isis = [] |
|---|
| 482 | | for idx,id in enumerate(self.id_list): |
|---|
| 483 | | isis.append(self.spiketrains[id].isi()) |
|---|
| 484 | | if not display: |
|---|
| 485 | | return isis |
|---|
| 486 | | else: |
|---|
| 487 | | ISI = numpy.array([]) |
|---|
| 488 | | for idx in xrange(self.N()): |
|---|
| 489 | | ISI = numpy.concatenate((ISI,isis[idx])) |
|---|
| 490 | | values, xaxis = numpy.histogram(ISI, nbins, normed=1) |
|---|
| 491 | | pylab.figure() |
|---|
| 492 | | pylab.plot(xaxis, values) |
|---|
| 493 | | pylab.xlabel("Inter Spike Interval (ms)", size="x-large") |
|---|
| 494 | | pylab.ylabel("% of Neurons", size="x-large") |
|---|
| 495 | | |
|---|
| 496 | | def isi_hist(self, bins): |
|---|
| 497 | | """Return the histogram of the ISI. |
|---|
| 498 | | |
|---|
| 499 | | bins may either be an integer, giving the number of bins (between the |
|---|
| 500 | | min and max of the data) or a list/array containing the lower edges of |
|---|
| 501 | | the bins. |
|---|
| 502 | | |
|---|
| 503 | | Returns a tuple, (histogram values, bin edges) |
|---|
| 504 | | """ |
|---|
| 505 | | isis = numpy.concatenate(self.isi()) |
|---|
| 506 | | return numpy.histogram(isis, bins=bins) |
|---|
| 507 | | |
|---|
| 508 | | |
|---|
| 509 | | def cv_isi(self, nbins=100, display=False): |
|---|
| 510 | | """ |
|---|
| 511 | | Return the list of all the cv coefficients for all the SpikeTrains objects |
|---|
| 512 | | within the SpikeList. If display is True, then it plots the distribution |
|---|
| 513 | | of the CVs |
|---|
| 514 | | """ |
|---|
| 515 | | cvs_isi = [] |
|---|
| 516 | | for idx,id in enumerate(self.id_list): |
|---|
| 517 | | isi = self.spiketrains[id].isi() |
|---|
| 518 | | if len(isi) > 1: |
|---|
| 519 | | cvs_isi.append(numpy.std(isi)/numpy.mean(isi)) |
|---|
| 520 | | if not display: |
|---|
| 521 | | return cvs_isi |
|---|
| 522 | | else: |
|---|
| 523 | | CV = numpy.array([]) |
|---|
| 524 | | for idx in xrange(len(cvs_isi)): |
|---|
| 525 | | CV = numpy.concatenate((CV,[cvs_isi[idx]])) |
|---|
| 526 | | values, xaxis = numpy.histogram(CV, nbins, normed=1) |
|---|
| 527 | | pylab.figure() |
|---|
| 528 | | pylab.plot(xaxis, values) |
|---|
| 529 | | pylab.xlabel("Inter Spike Interval CV", size="x-large") |
|---|
| 530 | | pylab.ylabel("% of Neurons", size="x-large") |
|---|
| 531 | | |
|---|
| 532 | | def cv_isi_hist(self, bins): |
|---|
| 533 | | cvs = numpy.array(self.cv_isi()) |
|---|
| 534 | | return numpy.histogram(cvs, bins=bins) |
|---|
| 535 | | |
|---|
| 536 | | def time_axis(self, time_bin): |
|---|
| 537 | | return numpy.arange(self.t_start, self.t_stop, time_bin) |
|---|
| 538 | | |
|---|
| 539 | | def mean_rate(self, t_start=None, t_stop=None): |
|---|
| 540 | | """ |
|---|
| 541 | | Return the mean firing rate averaged accross all SpikeTrains |
|---|
| 542 | | |
|---|
| 543 | | see mean_rates |
|---|
| 544 | | """ |
|---|
| 545 | | return numpy.mean(self.mean_rates(t_start, t_stop)) |
|---|
| 546 | | |
|---|
| 547 | | def mean_rate_std(self, t_start=None, t_stop=None): |
|---|
| 548 | | """ |
|---|
| 549 | | Std deviation of the Mean firing rate averaged accross all SpikeTrains |
|---|
| 550 | | |
|---|
| 551 | | see mean_rate |
|---|
| 552 | | """ |
|---|
| 553 | | return numpy.std(self.mean_rates(t_start, t_stop)) |
|---|
| 554 | | |
|---|
| 555 | | def mean_rates(self, t_start=None, t_stop=None): |
|---|
| 556 | | """ returns a vector of the size of id_list giving the mean rate for each neuron |
|---|
| 557 | | |
|---|
| 558 | | see SpikeTrain.mean_rate |
|---|
| 559 | | """ |
|---|
| 560 | | rates = [] |
|---|
| 561 | | for id in self.id_list: |
|---|
| 562 | | rates.append(self.spiketrains[id].mean_rate(t_start, t_stop)) |
|---|
| 563 | | |
|---|
| 564 | | return rates |
|---|
| 565 | | |
|---|
| 566 | | def rate_distribution(self, nbins=25, normalize=True, display=False): |
|---|
| 567 | | """ |
|---|
| 568 | | Return a vector with all the mean firing rates for all SpikeTrains. |
|---|
| 569 | | If display is True, then it plots the distribution of the rates |
|---|
| 570 | | """ |
|---|
| 571 | | #rates = numpy.zeros(self.N(), float) |
|---|
| 572 | | #for idx,id in enumerate(self.id_list): |
|---|
| 573 | | # rates[idx] = self.spiketrains[id].mean_firing_rate() |
|---|
| 574 | | rates = self.mean_rates() |
|---|
| 575 | | if not display: |
|---|
| 576 | | return rates |
|---|
| 577 | | else: |
|---|
| 578 | | values, xaxis = numpy.histogram(rates, nbins, normed=1) |
|---|
| 579 | | pylab.plot(xaxis,values) |
|---|
| 580 | | pylab.ylabel("% of Neurons", size="x-large") |
|---|
| 581 | | pylab.xlabel("Average Firing Rate (Hz)", size="x-large") |
|---|
| 582 | | |
|---|
| 583 | | def spike_histogram(self, time_bin, normalized=False, display=False): |
|---|
| 584 | | """ |
|---|
| 585 | | Generate an array with all the spike_histograms of all the SpikeTrains |
|---|
| 586 | | objects within the SpikeList. If display is True, then it plots the |
|---|
| 587 | | mean firing rate of the all population along time |
|---|
| 588 | | """ |
|---|
| 589 | | nbins = numpy.ceil((self.t_stop-self.t_start)/time_bin) |
|---|
| 590 | | firing_rate = numpy.zeros((self.N(),nbins), float) |
|---|
| 591 | | print "nbins = %d" % nbins |
|---|
| 592 | | for idx,id in enumerate(self.id_list): |
|---|
| 593 | | print idx, id |
|---|
| 594 | | firing_rate[idx,:] = self.spiketrains[id].time_histogram(time_bin,normalized) |
|---|
| 595 | | if not display: |
|---|
| 596 | | return firing_rate |
|---|
| 597 | | else: |
|---|
| 598 | | pylab.figure() |
|---|
| 599 | | pylab.plot(self.time_axis(time_bin),sum(firing_rate)/self.N()) |
|---|
| 600 | | pylab.ylabel("Mean Number of Spikes per bin", size="x-large") |
|---|
| 601 | | pylab.xlabel("Time (ms)", size="x-large") |
|---|
| 602 | | |
|---|
| 603 | | def firing_rate(self, time_bin, display=False): |
|---|
| 604 | | """ |
|---|
| 605 | | Calculate firing rate traces (in Hz) from arrays of spike times. |
|---|
| 606 | | |
|---|
| 607 | | Spike times and binwidth should are in milliseconds. |
|---|
| 608 | | Returns a tuple (list_of_firing_rate_traces,bins) |
|---|
| 609 | | |
|---|
| 610 | | >>> pylab.plot(output[0].time_axis(dt),sum(output.firing_rate(dt))) |
|---|
| 611 | | """ |
|---|
| 612 | | return self.spike_histogram(time_bin, normalized=True, display=display) |
|---|
| 613 | | |
|---|
| 614 | | # Compute the Fano Factor of the population. Need to be checked |
|---|
| 615 | | def fano_factor(self, time_bin): |
|---|
| 616 | | firing_rate = self.spike_histogram(time_bin) |
|---|
| 617 | | firing_rate = sum(firing_rate) |
|---|
| 618 | | fano = numpy.var(firing_rate)/numpy.mean(firing_rate) |
|---|
| 619 | | return fano |
|---|
| 620 | | |
|---|
| 621 | | def fano_factors_isi(self): |
|---|
| 622 | | """ returns a list containing the fano factors for each neuron""" |
|---|
| 623 | | fano_factors = [] |
|---|
| 624 | | for id in self.id_list: |
|---|
| 625 | | try: |
|---|
| 626 | | fano_factors.append(self.spiketrains[id].fano_factor_isi()) |
|---|
| 627 | | except: |
|---|
| 628 | | pass |
|---|
| 629 | | |
|---|
| 630 | | return fano_factors |
|---|
| 631 | | |
|---|
| 632 | | |
|---|
| 633 | | def id2position(self, id, dims): |
|---|
| 634 | | """ |
|---|
| 635 | | Allow to translate a gid (ranging from 0 to N*M) into |
|---|
| 636 | | a position on a N*M grid. dims should be given as a tuple |
|---|
| 637 | | (N,M) |
|---|
| 638 | | """ |
|---|
| 639 | | # Then we translate it in 2D |
|---|
| 640 | | if len(dims) == 1: |
|---|
| 641 | | return id |
|---|
| 642 | | if len(dims) == 2: |
|---|
| 643 | | x = numpy.floor(id/dims[0]) |
|---|
| 644 | | y = id % dims[1] |
|---|
| 645 | | return (x,y) |
|---|
| 646 | | |
|---|
| 647 | | |
|---|
| 648 | | def activity_map(self, dims, bounds=None, display=False): |
|---|
| 649 | | """ |
|---|
| 650 | | Generate a map of the activity during t_start and t_stop. |
|---|
| 651 | | If dims is a tuple, then cells are placed on a grid of size |
|---|
| 652 | | (N,M), else if dims is an array of size (2,nb_cells) with the |
|---|
| 653 | | x (first line) and y (second line) flotting positions of the cells, |
|---|
| 654 | | we generate a scatter plot. bounds is a parameters allowing to specify |
|---|
| 655 | | the range of the colorbar |
|---|
| 656 | | """ |
|---|
| 657 | | if isinstance(dims, tuple) or isinstance(dims, list): |
|---|
| 658 | | activity_map = numpy.zeros(dims,float) |
|---|
| 659 | | rates = self.mean_rates() |
|---|
| 660 | | for id in self.id_list: |
|---|
| 661 | | position = self.id2position(id, dims) |
|---|
| 662 | | activity_map[position] = rates[id] |
|---|
| 663 | | if not display: |
|---|
| 664 | | return activity_map |
|---|
| 665 | | else: |
|---|
| 666 | | pylab.figure() |
|---|
| 667 | | pylab.imshow(activity_map,interpolation='bicubic') |
|---|
| 668 | | pylab.colorbar() |
|---|
| 669 | | pylab.show() |
|---|
| 670 | | if bounds is not None: |
|---|
| 671 | | pylab.clim(bounds) |
|---|
| 672 | | elif isinstance(dims, numpy.ndarray): |
|---|
| 673 | | if not len(self.id_list) == len(dims[0]): |
|---|
| 674 | | raise Exception("Error, the number of positions does not match the number of cells in the SpikeList") |
|---|
| 675 | | rates = self.mean_rates() |
|---|
| 676 | | #for id in self.id_list: |
|---|
| 677 | | # rates.append(self.spiketrains[id].mean_firing_rate()) |
|---|
| 678 | | if not display: |
|---|
| 679 | | return rates |
|---|
| 680 | | else: |
|---|
| 681 | | x = dims[0,:] |
|---|
| 682 | | y = dims[1,:] |
|---|
| 683 | | pylab.scatter(x,y,c=rates) |
|---|
| 684 | | pylab.colorbar() |
|---|
| 685 | | if bounds is not None: |
|---|
| 686 | | pylab.clim(bounds) |
|---|
| 687 | | |
|---|
| 688 | | def pairwise_correlations(self, pairs, time_bin=1., display=False): |
|---|
| 689 | | """ |
|---|
| 690 | | Function to generate an array of cross correlation between computed |
|---|
| 691 | | between pairs of cells within the SpikeTrains. pairs should be therefore |
|---|
| 692 | | a list of (cell_id_1, cell_id_2). If display is True, then if plots the |
|---|
| 693 | | averaged cross correlation over all those pairs |
|---|
| 694 | | """ |
|---|
| 695 | | if len(pairs[0]) != len(pairs[1]): |
|---|
| 696 | | raise Exception("Pairs should have the same number of elements") |
|---|
| 697 | | nb_pairs = len(pairs[0]) |
|---|
| 698 | | spk1 = self.idsubSpikeList(pairs[0]) |
|---|
| 699 | | spk2 = self.idsubSpikeList(pairs[1]) |
|---|
| 700 | | hist_1 = spk1.spike_histogram(time_bin) |
|---|
| 701 | | hist_2 = spk2.spike_histogram(time_bin) |
|---|
| 702 | | print spk1, spk2 |
|---|
| 703 | | length = 2*len(hist_1[0]) |
|---|
| 704 | | results = numpy.zeros((nb_pairs,length), float) |
|---|
| 705 | | for idx in xrange(nb_pairs): |
|---|
| 706 | | # We need to avoid empty spike histogram, otherwise the ccf function |
|---|
| 707 | | # will give a nan vector |
|---|
| 708 | | if sum(hist_1[idx]) > 0 and sum(hist_2[idx] > 0): |
|---|
| 709 | | results[idx,:] = ccf(hist_1[idx],hist_2[idx]) |
|---|
| 710 | | if (display): |
|---|
| 711 | | results = sum(results)/nb_pairs |
|---|
| 712 | | pylab.figure() |
|---|
| 713 | | xaxis = time_bin*numpy.arange(-len(results)/2, len(results)/2) |
|---|
| 714 | | pylab.plot(xaxis, results) |
|---|
| 715 | | pylab.xlabel("Time (ms)", size="x-large") |
|---|
| 716 | | pylab.ylabel("Cross Correlation", size="x-large") |
|---|
| 717 | | else: |
|---|
| 718 | | return results |
|---|
| 719 | | |
|---|
| 720 | | def mean_rate_variance(self, time_bin): |
|---|
| 721 | | """ |
|---|
| 722 | | Function to extract the variance of the firing rate along time, |
|---|
| 723 | | if events are binned with a time bin. |
|---|
| 724 | | """ |
|---|
| 725 | | firing_rate = self.firing_rate(time_bin) |
|---|
| 726 | | return numpy.var(sum(firing_rate)/self.N()) |
|---|
| 727 | | |
|---|
| 728 | | def mean_rate_covariance(self, SpkList, time_bin): |
|---|
| 729 | | """ |
|---|
| 730 | | Function to extract the covariance of the firing rate along time, |
|---|
| 731 | | if events are binned with a time bin. We need a second SpikeList |
|---|
| 732 | | object with same time parameters to calculate the covariance |
|---|
| 733 | | """ |
|---|
| 734 | | if not isinstance(SpkList, SpikeList): |
|---|
| 735 | | raise Exception("Error, argument should be a SpikeList object") |
|---|
| 736 | | if not SpkList.get_time_parameters() == self.get_time_parameters(): |
|---|
| 737 | | raise Exception("Error, both SpikeLists should share common t_start, t_stop and dt") |
|---|
| 738 | | frate_1 = self.firing_rate(time_bin) |
|---|
| 739 | | frate_1 = sum(frate_1)/self.N() |
|---|
| 740 | | frate_2 = SpkList.firing_rate(time_bin) |
|---|
| 741 | | frate_2 = sum(frate_2)/SpkList.N() |
|---|
| 742 | | N = len(frate_1) |
|---|
| 743 | | cov = sum(frate_1*frate_2)/N-sum(frate_1)*sum(frate_2)/(N*N) |
|---|
| 744 | | return cov |
|---|
| 745 | | |
|---|
| 746 | | def raster_plot(self, id_list=None, t_start=None, t_stop=None, colors='b', subplot=None, size=1): |
|---|
| 747 | | """ |
|---|
| 748 | | Generate a raster plot of a certain number of cells in the |
|---|
| 749 | | SpikeList object. If id_list is an integer, then N ids will be randomly choosen |
|---|
| 750 | | in id_list. If this is a list, those id will be selected. |
|---|
| 751 | | Raster is made between t_start and t_stop (region of interest, not the global ones) |
|---|
| 752 | | and colors can be a list of color (each for one cells) or a single string |
|---|
| 753 | | to apply the same color to all the raster plots. |
|---|
| 754 | | """ |
|---|
| 755 | | if id_list == None: |
|---|
| 756 | | id_list = self.id_list |
|---|
| 757 | | spk = self.idsubSpikeList(id_list) |
|---|
| 758 | | elif id_list != self.id_list: |
|---|
| 759 | | spk = self.idsubSpikeList(id_list) |
|---|
| 760 | | id_list = spk.id_list |
|---|
| 761 | | else: |
|---|
| 762 | | spk = self.idsubSpikeList(id_list) |
|---|
| 763 | | |
|---|
| 764 | | if t_start is None: t_start = self.t_start |
|---|
| 765 | | if t_stop is None: t_stop = self.t_stop |
|---|
| 766 | | |
|---|
| 767 | | if subplot is None: |
|---|
| 768 | | pylab.figure() |
|---|
| 769 | | subplot = pylab |
|---|
| 770 | | if isinstance(colors, str): # this is much, much faster than a different colour per row |
|---|
| 771 | | ids, spike_times = self.as_ids_times() |
|---|
| 772 | | idx = numpy.where((spike_times >= t_start) & (spike_times <= t_stop))[0] |
|---|
| 773 | | spike_times = spike_times[idx] |
|---|
| 774 | | ids = ids[idx] |
|---|
| 775 | | if len(spike_times) > 0: |
|---|
| 776 | | print "Plotting %d points for %s" % (len(spike_times), self.label) |
|---|
| 777 | | subplot.scatter(spike_times, ids, s=size, c=colors) |
|---|
| 778 | | elif len(colors) != len(id_list): |
|---|
| 779 | | raise Exception("The numbers of colors does not match the number of cells!") |
|---|
| 780 | | else: # this is very slow in interactive mode |
|---|
| 781 | | for count,id in enumerate(spk.id_list): |
|---|
| 782 | | st = spk.spiketrains[id] |
|---|
| 783 | | idx = numpy.where((st.spike_times >= t_start) & (st.spike_times <= t_stop))[0] |
|---|
| 784 | | spikes = st.spike_times[idx] |
|---|
| 785 | | if len(spikes) > 0: |
|---|
| 786 | | subplot.scatter(spikes,count*numpy.ones(len(spikes)), s=size, c=colors[count]) |
|---|
| 787 | | xlabel = "Time (ms)" |
|---|
| 788 | | ylabel = "Neuron #" |
|---|
| 789 | | if hasattr(subplot, 'xlabel'): |
|---|
| 790 | | subplot.xlabel(xlabel, size="x-large") |
|---|
| 791 | | subplot.ylabel(ylabel, size="x-large") |
|---|
| 792 | | else: |
|---|
| 793 | | subplot.set_xlabel(xlabel, size="x-large") |
|---|
| 794 | | subplot.set_ylabel(ylabel, size="x-large") |
|---|
| 795 | | subplot.axis([t_start-10, t_stop+10, -1, len(self)+1]) |
|---|
| 796 | | |
|---|
| 797 | | # Clearly not optimal, but at least this is working... The best way to |
|---|
| 798 | | # do that function would be to go through the sorted list of all the spike times |
|---|
| 799 | | def activity_movie(self, dims, time_bin=10, bounds = (0,20), output="animation.mpg", fps=10): |
|---|
| 800 | | time = self.t_start |
|---|
| 801 | | files = [] |
|---|
| 802 | | while (time < self.t_stop): |
|---|
| 803 | | subSpkList = self.timesubSpikeList(time, time + time_bin) |
|---|
| 804 | | activity_map = subSpkList.activity_map(dims, bounds, display=True) |
|---|
| 805 | | caption = time+time_bin/2 |
|---|
| 806 | | pylab.title("Time = %g ms" %caption) |
|---|
| 807 | | fname = "_tmp_%g.jpg" %caption |
|---|
| 808 | | print " Saving Frame", fname |
|---|
| 809 | | pylab.savefig(fname) |
|---|
| 810 | | pylab.close() |
|---|
| 811 | | files.append(fname) |
|---|
| 812 | | time += time_bin |
|---|
| 813 | | print 'Making movie %s - this make take a while' %output |
|---|
| 814 | | command = "mencoder 'mf://_tmp_*.jpg' -mf type=jpg:fps=%d -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o %s" %(fps,output) |
|---|
| 815 | | print command |
|---|
| 816 | | os.system(command) |
|---|
| 817 | | ## cleanup |
|---|
| 818 | | print "Clean up...." |
|---|
| 819 | | for fname in files: os.remove(fname) |
|---|
| 820 | | |
|---|
| 821 | | def pw_corr_pearson(self,edge,bins,number_of_neuron_pairs): |
|---|
| 822 | | """ |
|---|
| 823 | | TODO: document/test this function |
|---|
| 824 | | |
|---|
| 825 | | """ |
|---|
| 826 | | #bins = edge[1]-edge[0] |
|---|
| 827 | | cor = numpy.zeros((number_of_neuron_pairs,)) |
|---|
| 828 | | [neuron_ids_array, spike_times_array] = self.as_list_id_list_time() |
|---|
| 829 | | |
|---|
| 830 | | # Pairwise correlation |
|---|
| 831 | | neuron_ids_unique = numpy.unique(neuron_ids_array) |
|---|
| 832 | | |
|---|
| 833 | | if len(neuron_ids_unique)==0: |
|---|
| 834 | | return (0,0) |
|---|
| 835 | | |
|---|
| 836 | | for count in range(number_of_neuron_pairs): |
|---|
| 837 | | # draw two neuron radomly out of the neuron_id_unique pool |
|---|
| 838 | | neuron_1_tmp = neuron_ids_unique[numpy.floor(numpy.random.uniform()*len(neuron_ids_unique))] |
|---|
| 839 | | neuron_2_tmp = neuron_ids_unique[numpy.floor(numpy.random.uniform()*len(neuron_ids_unique))] |
|---|
| 840 | | # get the spike_times of the neurons |
|---|
| 841 | | spike_times_1_tmp = self.spiketrains[neuron_1_tmp].spike_times |
|---|
| 842 | | spike_times_2_tmp = self.spiketrains[neuron_2_tmp].spike_times |
|---|
| 843 | | |
|---|
| 844 | | # hist of the spiketrains |
|---|
| 845 | | n1_hist = numpy.histogram(spike_times_1_tmp,bins=bins,range=edge) |
|---|
| 846 | | n2_hist = numpy.histogram(spike_times_2_tmp,bins=bins,range=edge) |
|---|
| 847 | | |
|---|
| 848 | | # correlation |
|---|
| 849 | | # TODO: normalize the cor, look in 1.6 in Kumar et al. |
|---|
| 850 | | # bruederle: the function corrcoeff actually implements the definition in Kumar 1.6 |
|---|
| 851 | | cov = numpy.corrcoef(n1_hist[0],n2_hist[0])[1][0] |
|---|
| 852 | | #print n1_hist[0] |
|---|
| 853 | | #print n2_hist[0] |
|---|
| 854 | | |
|---|
| 855 | | # bruederle: the expression 'cov' is already, per definition, the pearson correlation coefficient |
|---|
| 856 | | # see http://en.wikipedia.org/wiki/Correlation#The_sample_correlation |
|---|
| 857 | | cor[count] = cov |
|---|
| 858 | | # these two versions have been existing here before |
|---|
| 859 | | #cor[count] = cov/numpy.sqrt(n1_hist[0].var()*n2_hist[0].var()) |
|---|
| 860 | | #cor = numpy.append(cor,numpy.corrcoef(n1_hist[0],n2_hist[0])[1][0]) |
|---|
| 861 | | |
|---|
| 862 | | cor_coef_mean = cor.mean() |
|---|
| 863 | | cor_coef_std = cor.std() |
|---|
| 864 | | return cor_coef_mean, cor_coef_std |
|---|
| 865 | | |
|---|
| 866 | | # methods for different output formats |
|---|
| 867 | | def as_ids_times(self, relative=False, quantized=False): |
|---|
| 868 | | """Returns (array of ids, array of times).""" |
|---|
| 869 | | spikes = numpy.concatenate([st.spike_times for st in self.spiketrains.itervalues()]) |
|---|
| 870 | | ids = numpy.concatenate([id*numpy.ones((len(st.spike_times),)) for id,st in self.spiketrains.iteritems()]) |
|---|
| 871 | | assert len(ids) == len(spikes) |
|---|
| 872 | | return ids, spikes |
|---|
| 873 | | |
|---|
| 874 | | def as_list_id_time(self, relative=False, quantized=False): |
|---|
| 875 | | """ |
|---|
| 876 | | Returns a list/tuple of (id,t) tuples (id in range(0,N)). |
|---|
| 877 | | |
|---|
| 878 | | """ |
|---|
| 879 | | spikes =[] |
|---|
| 880 | | for id in self.id_list: |
|---|
| 881 | | for time in self.spiketrains[id].spike_times : |
|---|
| 882 | | spikes.append((id, time)) |
|---|
| 883 | | return spikes |
|---|
| 884 | | |
|---|
| 885 | | def as_list_id_list_time(self, relative=False, quantized=False): |
|---|
| 886 | | """ |
|---|
| 887 | | Returns (list of ids, list of times). |
|---|
| 888 | | |
|---|
| 889 | | """ |
|---|
| 890 | | spikes = self.as_list_id_time() |
|---|
| 891 | | ids, times = [], [] |
|---|
| 892 | | for spike in spikes: |
|---|
| 893 | | ids.append(spike[0]) |
|---|
| 894 | | times.append(spike[1]) |
|---|
| 895 | | |
|---|
| 896 | | return [ids, times] |
|---|
| 897 | | |
|---|
| 898 | | def as_id_list_times(self, relative=False, quantized=False): |
|---|
| 899 | | pass |
|---|
| 900 | | |
|---|
| 901 | | def as_time_list_ids(self, relative=False, quantized=False): |
|---|
| 902 | | pass |
|---|
| 903 | | |
|---|
| 904 | | def as_2byN_array(self, relative=False, quantized=False): |
|---|
| 905 | | pass |
|---|
| 906 | | |
|---|
| 907 | | def as_spikematrix(self): |
|---|
| 908 | | """ Returns a list giving for each neuron how many spikes per |
|---|
| 909 | | time bin, usually zero everywhere, one if a spike. |
|---|
| 910 | | Useful for plots etc. |
|---|
| 911 | | |
|---|
| 912 | | """ |
|---|
| 913 | | return self.firing_rate(self.dt) |
|---|
| 914 | | |
|---|
| 915 | | def as_pyNN_SpikeArray(self): |
|---|
| 916 | | """ |
|---|
| 917 | | to use in conjonction with SpikeSourceArray neurons |
|---|
| 918 | | |
|---|
| 919 | | >>> pop = Population((10,),SpikeSourceArray, {'spike_times': sl.as_pyNN_SpikeArray() }) |
|---|
| 920 | | |
|---|
| 921 | | """ |
|---|
| 922 | | spike_array = list([ [] for i in range(self.N())]) |
|---|
| 923 | | for spike in spikes: |
|---|
| 924 | | spike_array[spike[0]].append(spike[1]) |
|---|
| 925 | | |
|---|
| 926 | | return spike_array |
|---|
| 927 | | |
|---|
| 928 | | """ |
|---|
| 929 | | SpikeLists functions |
|---|
| 930 | | |
|---|
| 931 | | """ |
|---|
| 932 | | |
|---|
| 933 | | def population2spikelist(output, N, dt , t_start, t_stop): |
|---|
| 934 | | """ |
|---|
| 935 | | TODO sl2pynn |
|---|
| 936 | | """ |
|---|
| 937 | | import os, tempfile |
|---|
| 938 | | |
|---|
| 939 | | tmpdir = tempfile.mkdtemp() |
|---|
| 940 | | output_filename=os.path.join(tmpdir,'output.gdf') |
|---|
| 941 | | output.printSpikes(output_filename)# |
|---|
| 942 | | output_DATA = loadSpikeList(output_filename, N= N, dt = dt, t_start=t_start, t_stop=t_stop) |
|---|
| 943 | | os.remove(output_filename) |
|---|
| 944 | | os.rmdir(tmpdir) |
|---|
| 945 | | return output_DATA |
|---|
| 946 | | |
|---|
| 947 | | def readFile(filename, sepchar = "\t", skipchar = '#'): |
|---|
| 948 | | """ |
|---|
| 949 | | Function to read data. Should be optimize to deal with errors in the file |
|---|
| 950 | | that may append sometimes, when nest1 does not save the entire last line. |
|---|
| 951 | | It assumes that the data have been produced by pyNN under the format |
|---|
| 952 | | time (ms) gids (for raster) |
|---|
| 953 | | value (unit) gids (for continuous recordings like Vm, current, conductances) |
|---|
| 954 | | """ |
|---|
| 955 | | myfile = open(filename, "r", 10000) |
|---|
| 956 | | contents = myfile.readlines() |
|---|
| 957 | | myfile.close() |
|---|
| 958 | | data = [] |
|---|
| 959 | | for line in contents: |
|---|
| 960 | | stripped_line = line.lstrip() |
|---|
| 961 | | if (len(stripped_line) != 0): |
|---|
| 962 | | if (stripped_line[0] != skipchar): |
|---|
| 963 | | items = stripped_line.split(sepchar) |
|---|
| 964 | | data.append([int(float(items[1])), float(items[0])]) |
|---|
| 965 | | return(data) |
|---|
| 966 | | |
|---|
| 967 | | def get_header(filename): |
|---|
| 968 | | cmd = '' |
|---|
| 969 | | f = open(filename, 'r') |
|---|
| 970 | | for line in f.readlines(): |
|---|
| 971 | | if line[0] == '#': |
|---|
| 972 | | cmd += line[1:].strip() + ';' |
|---|
| 973 | | f.close() |
|---|
| 974 | | header = {} |
|---|
| 975 | | exec cmd in None, header |
|---|
| 976 | | return header |
|---|
| 977 | | |
|---|
| 978 | | def loadSpikeList(filename, id_list=None, dt = None, t_start=None, t_stop=None): |
|---|
| 979 | | """ |
|---|
| 980 | | Returns a SpikeList object from the tmp file saved by PyNN. |
|---|
| 981 | | |
|---|
| 982 | | The input is the filename where pyNN stored the spike list and the discretization step dt. |
|---|
| 983 | | it returns the spike list represented as a list of events. Events are tuples (relative time |
|---|
| 984 | | since last event, neuron_id); neuron_id and time are integers |
|---|
| 985 | | |
|---|
| 986 | | All times in milliseconds. |
|---|
| 987 | | |
|---|
| 988 | | The PyNN format (with compatible_output=True) is: |
|---|
| 989 | | * comment lines containing header information, |
|---|
| 990 | | * then one line per event: absolute time in ms, GID |
|---|
| 991 | | (see function printSpikes in nest1.py) |
|---|
| 992 | | |
|---|
| 993 | | """ |
|---|
| 994 | | if id_list is None: # try to obtain id_list from file header |
|---|
| 995 | | header = get_header(filename) |
|---|
| 996 | | if 'first_id' in header: |
|---|
| 997 | | n_cells = int(header['last_id']) - int(header['first_id']) + 1 |
|---|
| 998 | | id_list = n_cells |
|---|
| 999 | | else: |
|---|
| 1000 | | raise Exception("id_list not provided and cannot be inferred from file header.") |
|---|
| 1001 | | if isinstance(id_list,int): # allows to just specify the number of neurons |
|---|
| 1002 | | id_list = range(id_list) |
|---|
| 1003 | | spikes = readFile(filename) |
|---|
| 1004 | | print "Read %d spikes from %s" % (len(spikes), filename) |
|---|
| 1005 | | return SpikeList(spikes, id_list, dt, t_start, t_stop, label=filename) |
|---|
| 1006 | | |
|---|
| 1007 | | |
|---|
| 1008 | | def loadConductanceTraceList(filename, id_list, dt = None, t_start=None, t_stop=None): |
|---|
| 1009 | | if isinstance(id_list,int): # allows to just specify the number of neurons |
|---|
| 1010 | | id_list = range(id_list) |
|---|
| 1011 | | analog_signals = readFile(filename) |
|---|
| 1012 | | return ConductanceTraceList(analog_signals, id_list, dt, t_start, t_stop) |
|---|
| 1013 | | |
|---|
| 1014 | | |
|---|
| 1015 | | def loadMembraneTraceList(filename, id_list, dt = None, t_start=None, t_stop=None): |
|---|
| 1016 | | if isinstance(id_list,int): # allows to just specify the number of neurons |
|---|
| 1017 | | id_list = range(id_list) |
|---|
| 1018 | | analog_signals = readFile(filename) |
|---|
| 1019 | | return MembraneTraceList(analog_signals, id_list, dt, t_start, t_stop) |
|---|
| 1020 | | |
|---|
| 1021 | | |
|---|
| 1022 | | def loadCurrentTraceList(filename, id_list, dt = None, t_start=None, t_stop=None): |
|---|
| 1023 | | if isinstance(id_list,int): # allows to just specify the number of neurons |
|---|
| 1024 | | id_list = range(id_list) |
|---|
| 1025 | | analog_signals = readFile(filename) |
|---|
