bmtk icon indicating copy to clipboard operation
bmtk copied to clipboard

Problem with spike trains in tutorials 2 and 3.

Open aaberbach opened this issue 4 years ago • 8 comments

When the PoissonSpikeGenerator is used in both tutorials, the numbers passed in for times and firing_rate are in the wrong units. The problem is that when the spikes are eventually played into the vecstim in the virtual cells, the times are interpreted as ms instead of seconds. Therefore, the start and stop that is passed into times for psg should be in ms and the firing_rate passed into the psg should be in kHz instead of Hz.

For example, in tutorial 2, the spike train is supposed to be 3 seconds long with spikes firing at 10 Hz. In the tutorial they use firing_rate=10.0 and times=(0.0, 3.0), but what they really want is firing_rate=0.01 and times=(0.0, 3000.0).

aaberbach avatar Jun 11 '20 15:06 aaberbach

Looks alright to me?

https://github.com/AllenInstitute/bmtk/blob/845feda4fe421363d3f24e86806cda1e716943af/bmtk/utils/reports/spike_trains/spike_trains.py#L151-L170

PoissonSpikeGenerator seems to take care of the conversion just fine. Does the output in the h5 file look wrong?

tjbanks avatar Jun 11 '20 16:06 tjbanks

The problem is that the spike times are confined within the tstart and tstop. Therefore, given that spike times are interpreted as ms in the virtual cell instead of as seconds, there will only be spikes between 0.0 ms and 3.0 ms instead of between 0.0 s and 3.0 s.

aaberbach avatar Jun 11 '20 17:06 aaberbach

Do you have a way of showing this? Everything appears correct on my end - input seconds into spike generator - output spikes in ms in h5 file.

tjbanks avatar Jun 11 '20 17:06 tjbanks

If you go to tutorial 2 where the spikes are generated, the dataframe of spike times is printed out. You can see there that all of the times are between 0 and 3.

aaberbach avatar Jun 11 '20 18:06 aaberbach

All is correct. seconds are provided in the notebook and ms are in the output file.

from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator

psg = PoissonSpikeGenerator(population='mthalamus')
psg.add(node_ids=range(10),  # Have 10 nodes to match mthalamus
        firing_rate=10.0,    # 10 Hz, we can also pass in a nonhomoegenous function/array
        times=(0.0, 3.0))    # Firing starts at 0 s up to 3 s
psg.to_sonata('sim_ch02/inputs/mthalamus_spikes.h5')
(base) PS C:\Users\Tyler\Desktop\git_stage\bmtk\docs\tutorial\sim_ch02\inputs> python
Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)] :: Anaconda, Inc. on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import h5py
>>> h = h5py.File('mthalamus_spikes.h5')
>>> h['spikes']['mthalamus']['timestamps'].value
array([  30.66750922,  269.11768975,  408.63443826,  434.50709074,
        504.54228925,  558.7706086 ,  592.30142255,  912.12970708,
       1049.28120522, 1119.54991702, 1342.81833934, 1468.75187376,
       1474.2661857 , 1489.1999796 , 1492.49756877, 1551.39426571,
       1610.78904545, 1632.45015293, 1703.02358577, 2101.66411639,
       2198.45376484, 2355.29897342, 2530.37518602, 2541.69950934,
       2665.12894438, 2724.41874513, 2763.65720292, 2830.10398158,
       2940.67640501,   54.81065968,  258.4683457 ,  296.84502644,
        400.60894048,  559.97302337,  718.8627887 ,  720.79037251,
        784.32846776, 1112.47917087, 1196.8164137 , 1220.20566854,
       1325.73155772, 1381.9506117 , 1425.82385964, 1561.72811512,
       1658.36091854, 1823.01757338, 1906.29826343, 2016.63842227,
       2109.25472666, 2268.35522041, 2412.55512374, 2442.99586936,
       2490.45981633, 2604.81603513, 2620.23567554, 2716.22397924,
       2742.68135266, 2776.33468339, 2929.1422949 , 2946.69690089,
         25.45443301,   78.29100026,  210.58239176,  220.94374648,
        352.62848594,  715.8903775 ,  745.96706821,  967.54121267,
       1034.49421868, 1421.02064844, 1466.2484627 , 1530.75379061,
       1601.56413809, 1868.19851742, 1965.03626349, 2007.25755036,
       2090.7825628 , 2287.13976444, 2291.24810629, 2324.68644875,
       2427.29975022, 2535.48951982, 2542.02644274, 2577.54481858,
       2619.11000521, 2717.40150371, 2757.97128302, 2845.14981464,
       2881.14225588,   72.06387183,  292.46841714,  301.28184334,
        353.23273634,  366.07567202,  383.08321741,  630.67097046,
        644.19057521,  667.14083414,  943.42243813, 1094.89436121,
       1311.44310158, 1481.25556338, 1533.9342706 , 1767.65185973,
       1806.26412454, 1809.08265882, 1948.46959002, 1979.30657791,
       2030.92649472, 2062.74065821, 2065.508331  , 2231.23155166,
       2231.44402044, 2236.9710166 , 2353.97166587, 2354.15450939,
       2502.51170473, 2890.99748264, 2921.58112364,  122.38637209,
        160.43377781,  326.73395318,  382.20973969,  386.9266155 ,
        481.5329043 ,  507.33178744,  516.49123216,  535.3967916 ,
        576.8864983 ,  633.30455189,  699.00824779,  879.73323586,
        914.33594501,  942.67745391, 1140.7093397 , 1157.19050779,
       1161.25901751, 1353.19543387, 1356.10578727, 1436.52281621,
       1493.59122934, 1531.96952169, 1921.16308753, 2110.31417357,
       2149.81154741, 2159.36429698, 2161.21558378, 2415.51044685,
       2433.56103708, 2485.90495487, 2509.83459371,  234.78928567,
        437.90005729,  522.49537575,  567.80623141,  604.11242546,
        623.81190962,  696.65770259,  744.76765373, 1277.96598929,
       1604.49614738, 1732.94310478, 1809.8924921 , 1816.25057467,
       1857.7228066 , 2090.30842755, 2108.72551106, 2116.10981012,
       2174.44271532, 2266.54297279, 2428.30284799, 2749.22286306,
       2760.35274008, 2764.598118  , 2856.66934221, 2865.6542397 ,
        106.7851983 ,  109.31081292,  126.28869438,  181.92860133,
        197.30368934,  216.84708642,  321.67301859,  344.60527955,
        483.00650232,  521.82939177,  589.11221992,  692.65112378,
        835.11316597,  909.64377773,  916.67233529, 1060.11620204,
       1099.24424224, 1265.58281487, 1558.87232326, 1591.58273904,
       1605.1847068 , 1662.97958302, 1663.82007497, 1670.77361002,
       1737.60118458, 1791.76733784, 1870.24384606, 1959.89296596,
       2015.33082182, 2154.59137796, 2212.07042817, 2239.68542643,
       2717.411851  , 2801.11489733, 2848.0979937 , 2983.84342478,
       2985.35784082,   33.34773632,  112.40089013,  136.37592985,
        287.19334566,  314.4652381 ,  351.62114411,  499.93662804,
        563.8404097 ,  641.28962851,  643.33050428,  670.99820003,
        796.3579204 ,  850.13474908, 1268.90119526, 1340.2011028 ,
       1367.94100514, 1722.12395397, 1758.12046359, 1809.63470097,
       1814.82455382, 2084.21217228, 2091.19917278, 2199.30481455,
       2218.98361376, 2373.14825382, 2385.87983429, 2453.30359061,
       2554.7261475 , 2608.04525167, 2655.71634155, 2676.22941732,
       2694.93078168, 2911.56256151, 2924.68955943, 2935.14748157,
       2940.43061378, 2963.39977111,   61.91665604,  110.57369485,
        198.66546668,  206.69702413,  241.34697047,  353.49507223,
        379.75333069,  423.1497388 ,  722.81334078,  755.59589505,
        838.10288246,  991.4608314 , 1038.7287561 , 1139.60638735,
       1326.09379978, 1406.99722827, 1467.13626353, 1716.8118457 ,
       1864.26681844, 2000.84613904, 2005.91524541, 2247.7754892 ,
       2321.35547836, 2359.40556398, 2452.88826486, 2566.35008061,
       2725.77788655, 2725.96191775, 2923.92956009,   27.40333667,
        108.03882721,  221.30595966,  298.05707677,  336.81090206,
        489.42423283,  561.24225933,  564.25619141,  622.25489138,
        660.75738869,  664.12204272,  845.59767496,  912.93188608,
       1056.28089417, 1077.03687967, 1426.28996932, 1552.27784473,
       1568.63238409, 1679.18200673, 1693.54775199, 1727.09819347,
       1829.16535443, 1872.67621841, 2036.66438875, 2580.19393493,
       2623.07818558, 2705.20845362, 2856.37922314, 2939.56057952,
       2940.72698693, 2962.32777484, 2964.03140609, 2998.02433669])

tjbanks avatar Jun 11 '20 18:06 tjbanks

That is strange because I don't get that.

from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator

psg = PoissonSpikeGenerator(population='mthalamus')
psg.add(node_ids=range(10),  # Have 10 nodes to match mthalamus
        firing_rate=10,    # 10 Hz, we can also pass in a nonhomoegenous function/array
        times=(0.0, 3.0))    # Firing starts at 0 s up to 3 s
psg.to_sonata('sim_ch02/inputs/mthalamus_spikes.h5')

hf = h5py.File('sim_ch02/inputs/mthalamus_spikes.h5')
hf['spikes']['mthalamus']['timestamps'].value

array([1.16889791e-01, 1.52156000e-01, 2.23272272e-01, 3.75511015e-01, 4.26832916e-01, 7.62522046e-01, 8.16683391e-01, 8.85415550e-01, 9.72277102e-01, 9.88742626e-01, 1.02023457e+00, 1.10537383e+00, 1.11498020e+00, 1.18950151e+00, 1.31782793e+00, 1.35015988e+00, 1.48704437e+00, 1.73726397e+00, 1.74734116e+00, 1.84329248e+00, 2.01011934e+00, 2.02516988e+00, 2.19306926e+00, 2.24275437e+00, 2.50538734e+00, 2.59570441e+00, 2.64525064e+00, 2.65357898e+00, 2.67387561e+00, 2.67639655e+00, 2.94804606e+00, 1.84026755e-01, 2.68478232e-01, 2.79319978e-01, 2.93154293e-01, 2.98367279e-01, 3.90910575e-01, 4.64548209e-01, 5.00191522e-01, 5.67283500e-01, 6.62326001e-01, 6.91853635e-01, 9.40909222e-01, 9.46808762e-01, 9.55689245e-01, 1.06413318e+00, 1.24291399e+00, 1.26135864e+00, 1.46923090e+00, 1.52388251e+00, 1.52985857e+00, 1.55023123e+00, 1.58062899e+00, 1.61191682e+00, 1.63010333e+00, 1.79750290e+00, 1.94733335e+00, 1.95305198e+00, 1.99177683e+00, 2.21407699e+00, 2.37184766e+00, 2.51991519e+00, 2.59497707e+00, 2.59591234e+00, 2.69037912e+00, 2.73553377e+00, 2.75877534e+00, 2.81258734e+00, 2.96404007e+00, 2.99371185e+00, 1.12569562e-02, 1.41452126e-01, 1.43923504e-01, 4.60176803e-01, 4.79813988e-01, 5.37563153e-01, 6.51966939e-01, 8.14125871e-01, 9.90119193e-01, 1.01727954e+00, 1.06493157e+00, 1.15844271e+00, 1.19977841e+00, 1.48994265e+00, 1.57272417e+00, 1.58191877e+00, 1.61694642e+00, 1.68549588e+00, 1.81748485e+00, 1.89941016e+00, 1.98483349e+00, 2.26212478e+00, 2.35462655e+00, 2.64540405e+00, 2.72665640e+00, 2.88322819e+00, 2.98506857e+00, 9.33137459e-03, 8.05122454e-02, 1.07171585e-01, 2.07582944e-01, 3.54693519e-01, 3.86391310e-01, 4.22627402e-01, 5.18874798e-01, 6.84868308e-01, 6.94731436e-01, 7.40176962e-01, 7.83561971e-01, 8.57035994e-01, 9.06034017e-01, 1.06738512e+00, 1.08083873e+00, 1.08229260e+00, 1.14135470e+00, 1.28567319e+00, 1.30238372e+00, 1.34666618e+00, 1.44055160e+00, 1.48979887e+00, 1.64067580e+00, 1.70543984e+00, 1.71425089e+00, 1.72050839e+00, 1.84825289e+00, 1.85008288e+00, 1.92238916e+00, 2.38016389e+00, 2.47448024e+00, 2.69805256e+00, 2.84275413e+00, 2.81952373e-03, 7.27963102e-02, 2.62984889e-01, 3.23773127e-01, 5.27674235e-01, 5.30173796e-01, 5.56064026e-01, 6.30335415e-01, 7.51271555e-01, 7.53041223e-01, 9.06530337e-01, 9.81833153e-01, 1.14072267e+00, 1.14566245e+00, 1.24377914e+00, 1.26423686e+00, 1.34569388e+00, 1.47704697e+00, 1.59341286e+00, 1.71424732e+00, 1.73471159e+00, 1.79270140e+00, 1.83286261e+00, 1.97936799e+00, 2.03636665e+00, 2.18692925e+00, 2.32996553e+00, 2.39701327e+00, 2.45489166e+00, 2.62461229e+00, 2.68484941e+00, 2.85764403e+00, 2.90979609e+00, 4.84480485e-03, 5.76192545e-02, 2.06916246e-01, 2.24488802e-01, 4.07735337e-01, 6.09702076e-01, 6.72425683e-01, 6.99143048e-01, 7.41297897e-01, 8.26916844e-01, 9.10426659e-01, 1.01442229e+00, 1.04195214e+00, 1.08553764e+00, 1.16322923e+00, 1.17688688e+00, 1.19611708e+00, 1.28159943e+00, 1.69946667e+00, 1.74995796e+00, 1.75273346e+00, 1.80852615e+00, 1.86792197e+00, 1.88101217e+00, 2.13462828e+00, 2.34339363e+00, 2.61516528e+00, 2.68980100e+00, 2.69676887e+00, 2.74117955e+00, 2.74714914e+00, 2.84631328e+00, 2.97037442e+00, 3.87957671e-02, 1.70193806e-01, 1.80745748e-01, 1.95548845e-01, 4.13524525e-01, 5.56148155e-01, 5.91176210e-01, 6.73102788e-01, 6.96980202e-01, 8.20520821e-01, 9.94288863e-01, 1.05838884e+00, 1.11619482e+00, 1.21261886e+00, 1.37923636e+00, 1.42574991e+00, 1.46865565e+00, 1.51261199e+00, 1.53849238e+00, 1.65060840e+00, 1.65527293e+00, 1.67576370e+00, 1.71128426e+00, 1.73754292e+00, 1.81534682e+00, 1.86852959e+00, 1.95177068e+00, 2.08448197e+00, 2.10223045e+00, 2.26600276e+00, 2.40116574e+00, 2.45122925e+00, 2.50885718e+00, 2.59372234e+00, 2.62545415e+00, 2.68321708e+00, 2.80950680e+00, 2.88772044e+00, 2.92172647e+00, 1.06813619e-01, 1.12735072e-01, 1.85256492e-01, 3.12695024e-01, 6.12821556e-01, 6.59491929e-01, 7.13488759e-01, 7.39211042e-01, 8.43186076e-01, 9.37337183e-01, 1.10172776e+00, 1.22427139e+00, 1.30891131e+00, 1.40248383e+00, 1.51078028e+00, 1.53574006e+00, 1.53962982e+00, 1.61460961e+00, 1.74866325e+00, 1.76317903e+00, 1.77772475e+00, 1.84523523e+00, 1.96632875e+00, 2.00896571e+00, 2.07160900e+00, 2.10179924e+00, 2.12139191e+00, 2.20647755e+00, 2.22116080e+00, 2.52444609e+00, 2.93855602e+00, 1.18984617e-01, 1.43447967e-01, 2.02457232e-01, 2.49773571e-01, 2.67394507e-01, 3.03260610e-01, 3.70535232e-01, 3.90185025e-01, 8.20096476e-01, 8.44122446e-01, 1.02278805e+00, 1.02949491e+00, 1.27126857e+00, 1.43002265e+00, 1.72535772e+00, 1.72610225e+00, 1.79722320e+00, 1.84292441e+00, 1.87932084e+00, 1.96356298e+00, 2.03954551e+00, 2.20938665e+00, 2.29387472e+00, 2.37158973e+00, 2.50095359e+00, 2.53155639e+00, 2.53636417e+00, 2.67143815e+00, 2.67321999e+00, 2.71701099e+00, 2.75769338e+00, 2.83708953e+00, 2.89817489e+00, 2.98309510e+00, 1.71332037e-01, 5.09027304e-01, 5.44507576e-01, 6.04978159e-01, 6.07812449e-01, 7.48415169e-01, 8.63465953e-01, 8.85421857e-01, 9.11413458e-01, 1.15924629e+00, 1.20551859e+00, 1.25055729e+00, 1.41574895e+00, 1.59396283e+00, 1.63239518e+00, 1.85308442e+00, 2.13726165e+00, 2.24718757e+00, 2.31244453e+00, 2.43250355e+00, 2.49232628e+00, 2.52240785e+00, 2.77994554e+00, 2.81905002e+00, 2.94260566e+00])

aaberbach avatar Jun 11 '20 21:06 aaberbach

@tjbanks I get the same output as @aaberbach. I tested it both on NeuroVM (Ubuntu) and my own MacOS, same output.

from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator

psg = PoissonSpikeGenerator(population='mthalamus') psg.add(node_ids=range(10), # Have 10 nodes to match mthalamus firing_rate=10.0, # 10 Hz, we can also pass in a nonhomoegenous function/array times=(0.0, 3.0)) # Firing starts at 0 s up to 3 s

psg.to_sonata('sim_ch02/inputs/mthalamus_spikes.h5')

psg.to_dataframe()['timestamps'].max()

Out[9]: 2.9945507034122474

So all of the spikes get delivered in the first 3 ms of the simulation.

If I change it to

from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator

psg = PoissonSpikeGenerator(population='mthalamus') psg.add(node_ids=range(10), # Have 10 nodes to match mthalamus firing_rate=0.010, # 10 Hz, we can also pass in a nonhomoegenous function/array times=(0.0, 3000.0),) # Firing starts at 0 s up to 3 s

psg.to_sonata('sim_ch02/inputs/mthalamus_spikes.h5')

I get an output that makes sense. Could this be the issue: https://github.com/AllenInstitute/bmtk/blob/develop/bmtk/utils/reports/spike_trains/spike_trains.py#L130

latimerb avatar Jun 12 '20 21:06 latimerb

@latimerb @aaberbach - Got your output, my mistake. Was on an old commit.

It does look like you can supply units to the constructor psg = PoissonSpikeGenerator(population='mthalamus', units='ms')

but nothing is done with it here: https://github.com/AllenInstitute/bmtk/blob/develop/bmtk/utils/reports/spike_trains/spike_trains.py#L151-L169 or here: https://github.com/AllenInstitute/bmtk/blob/develop/bmtk/utils/reports/spike_trains/spike_train_buffer.py#L66

Perhaps there needs to be a *1000 somewhere if self.units == 'ms'?

It does look like quite a bit has changed since the commit overall... surely supplying the frequency as a float isn't expected behavior.

tjbanks avatar Jun 13 '20 03:06 tjbanks