## Plot spikes and network bursts


### Generate images from datasets and files from previous analyses containing spikes and network bursts.

__autoMEA__ can also be used to visualize activity from a dataset that was previously analyzed - for which spikes and bursts were already detected.

Here we will show how to, given a dataset containing signal, spikes and bursts, how to generate images to visualize the analyzed data, and compare the results with those obtained using __autoMEA__.


We first import __autoMEA__, together with __numpy__ and __pandas__, which will be used to manipulate data.

```python
import automea
import numpy as np 
import pandas as pd 
am = automea.Analysis()
```


We create a variable containing the path to the dataset, the dataset name, and the well we want to analyze in the dataset.

```python
path_to_dataset = '/Volumes/T7/automea_files_to_plot/'
dataset = 'MJ22C002-1782_DIV7.h5'
well = 'D5'
```


And we create separate folders to save signals and spikes rasters plots.

```python
import os  
plots_dir = f'plots/{dataset}_{well}'
if os.path.exists(plots_dir) is False:
    os.mkdir(plots_dir)
    os.mkdir(plots_dir+'/signals')
    os.mkdir(plots_dir+'/rasters')
```


We inform __autoMEA__ the path to find the dataset, and load it using the __loadh5()__ method.

```python
am.path_to_dataset = path_to_dataset
am.loadh5(dataset)
```


After loaded, we can check that __autoMEA__ is already storing a signal array.

```python
am.signal, am.signal.shape
```




    (array([[ 27,  23, -17, ..., -24,  31,  31],
            [-70,  15,  69, ...,  42,   4,  -8],
            [ -5, -25, -41, ..., -13, -54, -11],
            ...,
            [-17, -14,   5, ..., -17, -80,  17],
            [-36,   7,  -7, ..., -37, -19, -29],
            [  8, -20,  22, ...,   9,  73,  71]], dtype=int32),
     (237, 6000000))





Following, we use the __loadsignal()__ method, to keep in memory only the signal from the well we want to analyze.

```python
am.loadsignal(am.signal[am.wellsFromData == am.wellLabelIndexDict[well]])
```


Here we load a model, and use __autoMEA__ to detect all the variables we want to compare with the already-analyzed dataset.

```python
am.model_name = 'signal30.h5'
am.loadmodel()
```

    Metal device set to: Apple M1
    
    systemMemory: 8.00 GB
    maxCacheSize: 2.67 GB
    



```python
am.detect_threshold()
am.detect_spikes()
am.detect_reverbs(method = 'model')
am.detect_bursts()
am.detect_net('reverbs')
am.detect_net('bursts')
```


Now we load the previously-analyzed spikes and network bursts, from deticated 'csv' files. We filter the dataset to keep only the activity from the well we want to analyze, and save the spikes in an array named __spikes_manual__, and the network bursts in a __net_manual__ array.

```python
path_to_manual_files = '/Volumes/T7/automea_files_to_plot/Examples list/RHEB_DIV07/'
spikes_manual_file = 'Timestamps_MJ22C002_1782_DIV07_D5.csv'
netbursts_manual_file = 'NetworkBursts_MJ22C002_1782_DIV07_D5.csv'
```


```python
spikes_manual_df = pd.read_csv(path_to_manual_files + spikes_manual_file, encoding='unicode_escape')
```


```python
spikes_manual_df
```
<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Channel ID</th>
      <th>Channel Label</th>
      <th>Well ID</th>
      <th>Well Label</th>
      <th>Compound ID</th>
      <th>Compound Name</th>
      <th>Experiment</th>
      <th>Dose [pM]</th>
      <th>Dose Label</th>
      <th>Timestamp [µs]</th>
      <th>Maximum Amplitude [pV]</th>
      <th>Minimum Amplitude [pV]</th>
      <th>Peak-to-peak Amplitude [pV]</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>191</td>
      <td>21</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>426300</td>
      <td>16212288</td>
      <td>-25391304</td>
      <td>41603592</td>
    </tr>
    <tr>
      <th>1</th>
      <td>191</td>
      <td>21</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>435100</td>
      <td>5066340</td>
      <td>-11920800</td>
      <td>16987140</td>
    </tr>
    <tr>
      <th>2</th>
      <td>191</td>
      <td>21</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>452400</td>
      <td>7569708</td>
      <td>-14006940</td>
      <td>21576648</td>
    </tr>
    <tr>
      <th>3</th>
      <td>191</td>
      <td>21</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>494500</td>
      <td>7808124</td>
      <td>-12516840</td>
      <td>20324964</td>
    </tr>
    <tr>
      <th>4</th>
      <td>191</td>
      <td>21</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>503100</td>
      <td>50007756</td>
      <td>-59246376</td>
      <td>109254132</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>65010</th>
      <td>185</td>
      <td>34</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>598999200</td>
      <td>7629312</td>
      <td>-23185956</td>
      <td>30815268</td>
    </tr>
    <tr>
      <th>65011</th>
      <td>185</td>
      <td>34</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>599003300</td>
      <td>15079812</td>
      <td>-2264952</td>
      <td>17344764</td>
    </tr>
    <tr>
      <th>65012</th>
      <td>185</td>
      <td>34</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>599010600</td>
      <td>13470504</td>
      <td>-8702184</td>
      <td>22172688</td>
    </tr>
    <tr>
      <th>65013</th>
      <td>185</td>
      <td>34</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>599039800</td>
      <td>16391100</td>
      <td>-7092876</td>
      <td>23483976</td>
    </tr>
    <tr>
      <th>65014</th>
      <td>185</td>
      <td>34</td>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>599148400</td>
      <td>5721984</td>
      <td>-13887732</td>
      <td>19609716</td>
    </tr>
  </tbody>
</table>
<p>65015 rows × 13 columns</p>
</div>




```python
channel_labels = np.unique(spikes_manual_df[spikes_manual_df['Well Label'] == well]['Channel Label'].iloc[:].values)
channel_labels
```




    array([12, 13, 21, 22, 23, 24, 31, 32, 33, 34, 42, 43])




```python
spikes_manual = []
for channel_label in channel_labels:
    spikes_manual.append([])
    for index in spikes_manual_df[(spikes_manual_df['Channel Label'] == channel_label) & (spikes_manual_df['Well Label'] == well)].index:
        spikes_manual[-1].append(spikes_manual_df['Timestamp [µs]'][index]//100)
```


```python
netbursts_manual_df = pd.read_csv(path_to_manual_files + netbursts_manual_file, encoding='unicode_escape')
netbursts_manual_df

```
<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Well ID</th>
      <th>Well Label</th>
      <th>Compound ID</th>
      <th>Compound Name</th>
      <th>Experiment</th>
      <th>Dose [pM]</th>
      <th>Dose Label</th>
      <th>Start timestamp [µs]</th>
      <th>Duration [µs]</th>
      <th>Spike Count</th>
      <th>Spike Frequency [Hz]</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>494500</td>
      <td>250600</td>
      <td>296</td>
      <td>98.430</td>
    </tr>
    <tr>
      <th>1</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>14255400</td>
      <td>319800</td>
      <td>608</td>
      <td>158.432</td>
    </tr>
    <tr>
      <th>2</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>18646700</td>
      <td>444100</td>
      <td>749</td>
      <td>140.546</td>
    </tr>
    <tr>
      <th>3</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>26212300</td>
      <td>305800</td>
      <td>653</td>
      <td>177.949</td>
    </tr>
    <tr>
      <th>4</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>30838400</td>
      <td>265700</td>
      <td>603</td>
      <td>189.123</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>91</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>573895000</td>
      <td>311100</td>
      <td>216</td>
      <td>57.859</td>
    </tr>
    <tr>
      <th>92</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>578123900</td>
      <td>375300</td>
      <td>725</td>
      <td>160.982</td>
    </tr>
    <tr>
      <th>93</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>585406500</td>
      <td>363400</td>
      <td>487</td>
      <td>111.677</td>
    </tr>
    <tr>
      <th>94</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>591906400</td>
      <td>390700</td>
      <td>690</td>
      <td>147.172</td>
    </tr>
    <tr>
      <th>95</th>
      <td>22</td>
      <td>D5</td>
      <td>No compound</td>
      <td>No compound</td>
      <td>20181127_No_compound_m201118_DIV7_rescue4EB_001</td>
      <td>0</td>
      <td>Control</td>
      <td>598675200</td>
      <td>379000</td>
      <td>734</td>
      <td>161.390</td>
    </tr>
  </tbody>
</table>
<p>96 rows × 11 columns</p>
</div>




```python
net_manual = []
for index in netbursts_manual_df[netbursts_manual_df['Well Label']==well].index:
    start_ts = netbursts_manual_df['Start timestamp [µs]'][index]
    duration = netbursts_manual_df['Duration [µs]'][index]
    net_manual.append([start_ts//100, (start_ts+duration)//100])
```


After having spikes and network bursts both from a previously performed analysis, and obtained with __autoMEA__, we use the __plot_raster()__ method to visualize raster plots, so that the two version can be compared.

```python
duration = 60
duration_shorter = 5

for i, start_time in enumerate(np.random.randint(0,500,size=3)):

    end_time = start_time + duration

    am.plot_raster(am.spikes, net_bursts = am.net_reverbs, start_time = start_time, end_time = end_time, 
                   save = f'{plots_dir}/rasters/raster_start{start_time}_duration{duration}_allwell_prednetreverbs.eps', show = False)
    am.plot_raster(am.spikes, net_bursts = net_manual, start_time = start_time, end_time = end_time, 
                   save = f'{plots_dir}/rasters/raster_start{start_time}_duration{duration}_allwell_manual.eps', show = False)

    for j in range(5):
        start_time_shorter = start_time + np.random.randint(0,95)
        channel_random  = np.random.randint(len(am.signal))
        end_time_shorter = start_time_shorter + duration_shorter
        am.plot_raster(am.spikes[channel_random], net_bursts = am.net_reverbs, start_time = start_time_shorter, end_time = end_time_shorter, 
                        save = f'{plots_dir}/rasters/raster_start{start_time}_duration{duration}_zoomedstart{start_time+start_time_shorter}_zoomedduration{duration_shorter}_prednetreverbs.eps', show = False)
        am.plot_raster(am.spikes[channel_random], net_bursts = net_manual, start_time = start_time_shorter, end_time = end_time_shorter,
                        save = f'{plots_dir}/rasters/raster_start{start_time}_duration{duration}_zoomedstart{start_time+start_time_shorter}_zoomedduration{duration_shorter}_manual.eps', show = False)

```


And finally we plot a (randomly-selected) 5-second activity with signal and threshold, together with __autoMEA__ and manual network bursts, to compare them.

```python
duration_shorter = 5

for _ in range(20):
    random_net = am.net_reverbs[np.random.randint(len(am.net_reverbs))]
    middle_of_net = random_net[0] + ((random_net[-1] - random_net[0])/2)
    start_time = middle_of_net/am.samplingFreq - 2.5
    channel_random  = np.random.randint(len(am.signal))


    am.plot_window(am.signal[channel_random], start_time=start_time, duration=duration_shorter,
                    threshold=am.threshold[channel_random], spikes = am.spikes[channel_random],
                    bursts=am.net_reverbs, net_bursts=net_manual, yunits = 'V', show = False,
                    save = f'{plots_dir}/signals/signal_channel{channel_random}_start{start_time}_duration5.eps')

```