diff --git a/caiman/base/movies.py b/caiman/base/movies.py
index e94ddd218..e267abadb 100644
--- a/caiman/base/movies.py
+++ b/caiman/base/movies.py
@@ -13,11 +13,11 @@
 import ipywidgets as widgets
 import logging
 import matplotlib
+import matplotlib.pyplot as plt
 import numpy as np
 import os
 import pathlib
 import pims
-import pylab as pl
 import scipy
 import skimage
 import sklearn
@@ -893,8 +893,8 @@ def local_correlations(self,
                                             order_mean=order_mean)
                 Cn = np.maximum(Cn, rho)
                 if do_plot:
-                    pl.imshow(Cn, cmap='gray')
-                    pl.pause(.1)
+                    plt.imshow(Cn, cmap='gray')
+                    plt.pause(.1)
 
             logging.debug('number of chunks:' + str(n_chunks - 1) + ' frames: ' +
                           str([(n_chunks - 1) * frames_per_chunk, T]))
@@ -904,8 +904,8 @@ def local_correlations(self,
                                         order_mean=order_mean)
             Cn = np.maximum(Cn, rho)
             if do_plot:
-                pl.imshow(Cn, cmap='gray')
-                pl.pause(.1)
+                plt.imshow(Cn, cmap='gray')
+                plt.pause(.1)
 
         return Cn
 
@@ -1117,7 +1117,7 @@ def to_2D(self, order='F') -> np.ndarray:
         T = self.shape[0]
         return np.reshape(self, (T, -1), order=order)
 
-    def zproject(self, method: str = 'mean', cmap=pl.cm.gray, aspect='auto', **kwargs) -> np.ndarray:
+    def zproject(self, method: str = 'mean', cmap=matplotlib.cm.gray, aspect='auto', **kwargs) -> np.ndarray:
         """
         Compute and plot projection across time:
 
@@ -1141,7 +1141,7 @@ def zproject(self, method: str = 'mean', cmap=pl.cm.gray, aspect='auto', **kwarg
             zp = np.std(self, axis=0)
         else:
             raise Exception('Method not implemented')
-        pl.imshow(zp, cmap=cmap, aspect=aspect, **kwargs)
+        plt.imshow(zp, cmap=cmap, aspect=aspect, **kwargs)
         return zp
 
     def play(self,
@@ -1175,7 +1175,7 @@ def play(self,
             interpolation:
                 interpolation method for 'opencv' and 'embed_opencv' backends
 
-            backend: 'pylab', 'notebook', 'opencv' or 'embed_opencv'; the latter 2 are much faster
+            backend: 'opencv', 'embed_opencv', 'pyplot', 'notebook': the first two are much faster
 
             do_loop: Whether to loop the video
 
@@ -2363,7 +2363,7 @@ def play_movie(movie,
         interpolation:
             interpolation method for 'opencv' and 'embed_opencv' backends
 
-        backend: 'pylab', 'notebook', 'opencv' or 'embed_opencv'; the latter 2 are much faster
+        backend: 'opencv', 'embed_opencv', 'pyplot', or 'notebook': the first two are much faster
 
         do_loop: Whether to loop the video
 
@@ -2399,8 +2399,8 @@ def play_movie(movie,
     """
     # todo: todocument
     it = True if (isinstance(movie, list) or isinstance(movie, tuple) or isinstance(movie, str)) else False
-    if backend == 'pylab':
-        logging.warning('*** Using pylab. This might be slow. If you can use the opencv backend it may be faster')
+    if backend == 'pyplot':
+        logging.warning('Using pyplot back end: not recommended. Using opencv will yield faster, higher-quality results.')
 
     gain = float(gain)     # convert to float in case we were passed an int
     if q_max < 100:
@@ -2444,26 +2444,26 @@ def process_frame(iddxx, frame, bord_px, magnification, interpolation, minmov, m
                         thickness=1)
         return frame
 
-    if backend == 'pylab':
-        pl.ion()
-        fig = pl.figure(1)
+    if backend == 'pyplot':
+        plt.ion()
+        fig = plt.figure(1)
         ax = fig.add_subplot(111)
         ax.set_title("Play Movie")
         im = ax.imshow((offset + (load(movie, subindices=slice(0,2), var_name_hdf5=var_name_hdf5) if it else movie)[0] - minmov) * gain / (maxmov - minmov + offset),
-                       cmap=pl.cm.gray,
+                       cmap=plt.cm.gray,
                        vmin=0,
                        vmax=1,
-                       interpolation='none')                                            # Blank starting image
+                       interpolation='none')  # Blank starting image
         fig.show()
         im.axes.figure.canvas.draw()
-        pl.pause(1)
+        plt.pause(1)
 
     elif backend == 'notebook':
         # First set up the figure, the axis, and the plot element we want to animate
-        fig = pl.figure()
-        im = pl.imshow(next(load_iter(movie, subindices=slice(0,1), var_name_hdf5=var_name_hdf5))\
-                    if it else movie[0], interpolation='None', cmap=pl.cm.gray)
-        pl.axis('off')
+        fig = plt.figure()
+        im = plt.imshow(next(load_iter(movie, subindices=slice(0,1), var_name_hdf5=var_name_hdf5))\
+                    if it else movie[0], interpolation='None', cmap=matplotlib.cm.gray)
+        plt.axis('off')
 
         if it:
             m_iter = load_iter(movie, subindices, var_name_hdf5)
@@ -2508,7 +2508,7 @@ def view(button):
                     frame_sum = 0
                     display_handle.update(Image(data=cv2.imencode(
                             '.jpg', np.clip((frame * 255.), 0, 255).astype('u1'))[1].tobytes()))
-                    pl.pause(1. / fr)
+                    plt.pause(1. / fr)
                 if stopButton.value==True:
                     break
         display(stopButton)
@@ -2556,17 +2556,17 @@ def view(button):
                 elif backend == 'embed_opencv' and not save_movie:
                     break
 
-                elif backend == 'pylab':
+                elif backend == 'pyplot':
                     if bord_px is not None and np.sum(bord_px) > 0:
                         frame = frame[bord_px:-bord_px, bord_px:-bord_px]
                     im.set_data((offset + frame) * gain / maxmov)
                     ax.set_title(str(iddxx))
-                    pl.axis('off')
+                    plt.axis('off')
                     fig.canvas.draw()
-                    pl.pause(1. / fr * .5)
-                    ev = pl.waitforbuttonpress(1. / fr * .5)
+                    plt.pause(1. / fr * .5)
+                    ev = plt.waitforbuttonpress(1. / fr * .5)
                     if ev is not None:
-                        pl.close()
+                        plt.close()
                         break
 
                 elif backend == 'notebook':
diff --git a/caiman/base/timeseries.py b/caiman/base/timeseries.py
index d4d10cfb1..f8868cb1c 100644
--- a/caiman/base/timeseries.py
+++ b/caiman/base/timeseries.py
@@ -9,9 +9,9 @@
 from dateutil.tz import tzlocal
 import h5py
 import logging
+import matplotlib.pyplot as plt
 import numpy as np
 import os
-import pylab as plt
 from pynwb import NWBHDF5IO, NWBFile
 from pynwb.ophys import TwoPhotonSeries, OpticalChannel
 from pynwb.device import Device
diff --git a/caiman/base/traces.py b/caiman/base/traces.py
index 3e829e0ed..5e9962b79 100644
--- a/caiman/base/traces.py
+++ b/caiman/base/traces.py
@@ -3,8 +3,9 @@
 import cv2
 import logging
 import numpy as np
-import pylab as pl
-pl.ion()
+import matplotlib
+import matplotlib.pyplot as plt
+plt.ion()
 
 import caiman.base.timeseries
 
@@ -80,7 +81,7 @@ def computeDFF(self, window_sec=5, minQuantile=20):
     def resample(self, fx=1, fy=1, fz=1, interpolation=cv2.INTER_AREA):
         raise Exception('Not Implemented. Look at movie resize')
 
-    def plot(self, stacked=True, subtract_minimum=False, cmap=pl.cm.jet, **kwargs):
+    def plot(self, stacked=True, subtract_minimum=False, cmap=matplotlib.cm.jet, **kwargs):
         """Plot the data
 
         author: ben deverett
@@ -91,7 +92,7 @@ def plot(self, stacked=True, subtract_minimum=False, cmap=pl.cm.jet, **kwargs):
             subtract_minimum : bool
                 subtract minimum from each individual trace
             cmap : matplotlib.LinearSegmentedColormap
-                color map for display. Options are found in pl.colormaps(), and are accessed as pl.cm.my_favourite_map
+                color map for display. Options are found in matplotlib.colormaps(), and are accessed as matplotlib.cm.my_favourite_map
             kwargs : dict
                 any arguments accepted by matplotlib.plot
 
@@ -104,7 +105,7 @@ def plot(self, stacked=True, subtract_minimum=False, cmap=pl.cm.jet, **kwargs):
         if len(d.shape) > 1:
             n = d.shape[1]
 
-        ax = pl.gca()
+        ax = plt.gca()
 
         colors = cmap(np.linspace(0, 1, n))
         ax.set_color_cycle(colors)
@@ -122,7 +123,7 @@ def plot(self, stacked=True, subtract_minimum=False, cmap=pl.cm.jet, **kwargs):
         ax2.set_yticklabels([str(i) for i in range(n)], weight='bold')
         [l.set_color(c) for l, c in zip(ax2.get_yticklabels(), colors)]
 
-        pl.gcf().canvas.draw()
+        plt.gcf().canvas.draw()
         return ax
 
     def extract_epochs(self, trigs=None, tb=1, ta=1):
diff --git a/caiman/behavior/behavior.py b/caiman/behavior/behavior.py
index 6828c37d2..ce9fc87be 100644
--- a/caiman/behavior/behavior.py
+++ b/caiman/behavior/behavior.py
@@ -5,8 +5,9 @@
 """
 
 import cv2
+import matplotlib
+import matplotlib.pyplot as plt
 import numpy as np
-import pylab as pl
 import scipy
 from scipy.sparse import coo_matrix
 from scipy.io import loadmat
@@ -38,14 +39,14 @@ def select_roi(img: np.ndarray, n_rois: int = 1) -> list:
 
     masks = []
     for _ in range(n_rois):
-        fig = pl.figure()
-        pl.imshow(img, cmap=pl.cm.gray)
+        fig = plt.figure()
+        plt.imshow(img, cmap=matplotlib.cm.gray)
         pts = fig.ginput(0, timeout=0)
         mask = np.zeros(np.shape(img), dtype=np.int32)
         pts = np.asarray(pts, dtype=np.int32)
         cv2.fillConvexPoly(mask, pts, (1, 1, 1), lineType=cv2.LINE_AA)
         masks.append(mask)
-        pl.close()
+        plt.close()
 
     return masks
 
@@ -353,12 +354,12 @@ def extract_components(mov_tot,
 
 def plot_components(sp_filt, t_trace) -> None:
     # todo: todocument
-    pl.figure()
+    plt.figure()
     count = 0
     for comp, tr in zip(sp_filt, t_trace):
         count += 1
-        pl.subplot(6, 2, count)
-        pl.imshow(comp)
+        plt.subplot(6, 2, count)
+        plt.imshow(comp)
         count += 1
-        pl.subplot(6, 2, count)
-        pl.plot(tr)
+        plt.subplot(6, 2, count)
+        plt.plot(tr)
diff --git a/caiman/motion_correction.py b/caiman/motion_correction.py
index c6fd3dc28..a819fff95 100644
--- a/caiman/motion_correction.py
+++ b/caiman/motion_correction.py
@@ -41,11 +41,11 @@
 import gc
 import itertools
 import logging
+import matplotlib.pyplot as plt
 import numpy as np
 from numpy.fft import ifftshift
 import os
 import sys
-import pylab as pl
 import tifffile
 from typing import Optional
 from skimage.transform import resize as resize_sk
@@ -362,8 +362,8 @@ def motion_correct_pwrigid(self, save_movie:bool=True, template:np.ndarray=None,
                     indices=self.indices)
             if not self.is3D:
                 if show_template:
-                    pl.imshow(new_template_els)
-                    pl.pause(.5)
+                    plt.imshow(new_template_els)
+                    plt.pause(.5)  #TODO: figure out if pausing half-second is necessary
             if np.isnan(np.sum(new_template_els)):
                 raise Exception(
                     'Template contains NaNs, something went wrong. Reconsider the parameters')
@@ -877,10 +877,10 @@ def motion_correct_online(movie_iterable, add_to_movie, max_shift_w=25, max_shif
                     template = np.median(buffer_templates, 0)
 
                 if show_template:
-                    pl.cla()
-                    pl.imshow(template, cmap='gray', vmin=250,
+                    plt.cla()
+                    plt.imshow(template, cmap='gray', vmin=250,
                               vmax=350, interpolation='none')
-                    pl.pause(.001)
+                    plt.pause(.001)
 
                 logging.debug('Relative change in template:' + str(
                     np.sum(np.abs(template - template_old)) / np.sum(np.abs(template))))
@@ -2687,23 +2687,23 @@ def compute_metrics_motion_correction(fname, final_size_x, final_size_y, swap_di
                 cv2.putText(vid_frame, 'x_flow', (2 * dims[0] + 10, 20), fontFace=5, fontScale=0.8, color=(
                     0, 255, 0), thickness=1)
                 cv2.imshow('frame', vid_frame)
-                pl.pause(1 / fr_play)
+                cv2.waitKey(1 / fr_play) # to pause between frames
                 if cv2.waitKey(1) & 0xFF == ord('q'):
                     break
             else:
-                pl.subplot(1, 3, 1)
-                pl.cla()
-                pl.imshow(fr, vmin=m_min, vmax=m_max, cmap='gray')
-                pl.title('movie')
-                pl.subplot(1, 3, 3)
-                pl.cla()
-                pl.imshow(flow[:, :, 1], vmin=vmin, vmax=vmax)
-                pl.title('y_flow')
-                pl.subplot(1, 3, 2)
-                pl.cla()
-                pl.imshow(flow[:, :, 0], vmin=vmin, vmax=vmax)
-                pl.title('x_flow')
-                pl.pause(1 / fr_play)
+                plt.subplot(1, 3, 1)
+                plt.cla()
+                plt.imshow(fr, vmin=m_min, vmax=m_max, cmap='gray')
+                plt.title('movie')
+                plt.subplot(1, 3, 3)
+                plt.cla()
+                plt.imshow(flow[:, :, 1], vmin=vmin, vmax=vmax)
+                plt.title('y_flow')
+                plt.subplot(1, 3, 2)
+                plt.cla()
+                plt.imshow(flow[:, :, 0], vmin=vmin, vmax=vmax)
+                plt.title('x_flow')
+                plt.pause(1 / fr_play)
 
         n = np.linalg.norm(flow)
         flows.append(flow)
diff --git a/caiman/source_extraction/cnmf/utilities.py b/caiman/source_extraction/cnmf/utilities.py
index 532d7c271..c73678adb 100644
--- a/caiman/source_extraction/cnmf/utilities.py
+++ b/caiman/source_extraction/cnmf/utilities.py
@@ -13,7 +13,7 @@
 import numpy as np
 import os
 import pathlib
-import pylab as pl
+import matplotlib.pyplot as plt
 import scipy
 from scipy.sparse import spdiags, issparse, csc_matrix, csr_matrix
 import scipy.ndimage as ndi
@@ -749,7 +749,7 @@ def manually_refine_components(Y, xxx_todo_changeme, A, C, Cn, thr=0.9, display_
 
     x, y = np.mgrid[0:d1:1, 0:d2:1]
 
-    pl.imshow(Cn, interpolation=None, cmap=cmap)
+    plt.imshow(Cn, interpolation=None, cmap=cmap)
     cm = caiman.base.rois.com(A, d1, d2)
 
     Bmat = np.zeros((np.minimum(nr, max_number), d1, d2))
@@ -763,13 +763,13 @@ def manually_refine_components(Y, xxx_todo_changeme, A, C, Cn, thr=0.9, display_
 
     T = np.shape(Y)[-1]
 
-    pl.close()
-    fig = pl.figure()
-    ax = pl.gca()
+    plt.close()
+    fig = plt.figure()
+    ax = plt.gca()
     ax.imshow(Cn, interpolation=None, cmap=cmap,
               vmin=np.percentile(Cn[~np.isnan(Cn)], 1), vmax=np.percentile(Cn[~np.isnan(Cn)], 99))
     for i in range(np.minimum(nr, max_number)):
-        pl.contour(y, x, Bmat[i], [thr])
+        plt.contour(y, x, Bmat[i], [thr])
 
     if display_numbers:
         for i in range(np.minimum(nr, max_number)):
@@ -812,8 +812,8 @@ def manually_refine_components(Y, xxx_todo_changeme, A, C, Cn, thr=0.9, display_
             Bvec = np.zeros(d)
             Bvec[indx] = cumEn
             bmat = np.reshape(Bvec, np.shape(Cn), order='F')
-            pl.contour(y, x, bmat, [thr])
-            pl.pause(.01)
+            plt.contour(y, x, bmat, [thr])
+            plt.pause(.01)
 
         elif pts == []:
             break
diff --git a/caiman/tests/comparison_humans.py b/caiman/tests/comparison_humans.py
index e61664ff4..ea12f1fd4 100644
--- a/caiman/tests/comparison_humans.py
+++ b/caiman/tests/comparison_humans.py
@@ -23,10 +23,10 @@
     print('Not launched under iPython')
 
 import caiman as cm
+import matplotlib.pyplot as plt
 import numpy as np
 import os
 import time
-import pylab as pl
 import scipy
 import sys
 import logging
@@ -453,16 +453,14 @@
         'downsample_ratio': .2,
         'thr_plot': 0.8
     }
-
-    pl.rcParams['pdf.fonttype'] = 42
-    font = {'family': 'Arial',
-            'weight': 'regular',
-            'size': 20}
-    pl.rc('font', **font)
-
     plot_results = False
     if plot_results:
-        pl.figure(figsize=(30, 20))
+        plt.rcParams['pdf.fonttype'] = 42
+        font = {'family': 'Arial',
+                'weight': 'regular',
+                'size': 20}
+        plt.rc('font', **font)
+        plt.figure(figsize=(30, 20))
 
     tp_gt, tp_comp, fn_gt, fp_comp, performance_cons_off = compare_components(gt_estimate, cnm2.estimates,
                                                                               Cn=Cn_orig, thresh_cost=.8,
diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py
index 7bc2090ea..76d714655 100644
--- a/caiman/utils/visualization.py
+++ b/caiman/utils/visualization.py
@@ -11,10 +11,10 @@
 from math import sqrt, ceil
 import matplotlib
 import matplotlib.patches
+import matplotlib.pyplot as plt
 import matplotlib.widgets
 import numpy as np
 from numpy.typing import ArrayLike
-import pylab as pl
 from scipy.ndimage import center_of_mass, median_filter
 from scipy.sparse import issparse, spdiags, coo_matrix, csc_matrix
 from skimage.measure import find_contours
@@ -74,7 +74,7 @@ def view_patches(Yr, A, C, b, f, d1, d2, YrA=None, secs=1):
 
     """
 
-    pl.ion()
+    plt.ion()
     nr, T = C.shape
     nb = f.shape[0]
     A2 = A.copy()
@@ -88,37 +88,37 @@ def view_patches(Yr, A, C, b, f, d1, d2, YrA=None, secs=1):
 
     A = A.todense()
     bkgrnd = np.reshape(b, (d1, d2) + (nb,), order='F')
-    fig = pl.figure()
-    thismanager = pl.get_current_fig_manager()
+    fig = plt.figure()
+    thismanager = plt.get_current_fig_manager()
     thismanager.toolbar.pan()
     print('In order to scroll components you need to click on the plot')
     sys.stdout.flush()
     for i in range(nr + 1):
         if i < nr:
             ax1 = fig.add_subplot(2, 1, 1)
-            pl.imshow(np.reshape(np.array(A[:, i]) / nA2[i],
+            plt.imshow(np.reshape(np.array(A[:, i]) / nA2[i],
                                  (d1, d2), order='F'), interpolation='None')
             ax1.set_title('Spatial component ' + str(i + 1))
             ax2 = fig.add_subplot(2, 1, 2)
-            pl.plot(np.arange(T), np.squeeze(
+            plt.plot(np.arange(T), np.squeeze(
                 np.array(Y_r[i, :])), 'c', linewidth=3)
-            pl.plot(np.arange(T), np.squeeze(
+            plt.plot(np.arange(T), np.squeeze(
                 np.array(C[i, :])), 'r', linewidth=2)
             ax2.set_title('Temporal component ' + str(i + 1))
             ax2.legend(labels=['Filtered raw data', 'Inferred trace'])
 
             if secs > 0:
-                pl.pause(secs)
+                plt.pause(secs)
             else:
-                pl.waitforbuttonpress()
+                plt.waitforbuttonpress()
 
             fig.delaxes(ax2)
         else:
             ax1 = fig.add_subplot(2, 1, 1)
-            pl.imshow(bkgrnd[:, :, i - nr], interpolation='None')
+            plt.imshow(bkgrnd[:, :, i - nr], interpolation='None')
             ax1.set_title('Spatial background ' + str(i - nr + 1))
             ax2 = fig.add_subplot(2, 1, 2)
-            pl.plot(np.arange(T), np.squeeze(np.array(f[i - nr, :])))
+            plt.plot(np.arange(T), np.squeeze(np.array(f[i - nr, :])))
             ax2.set_title('Temporal background ' + str(i - nr + 1))
 
 
@@ -559,7 +559,7 @@ def nb_view_patches3d(Y_r, A, C, dims, image_type='mean', Yr=None,
         coors = [get_contours(proj_[i], tmp[i].shape, thr=thr)
                  for i in range(3)]
 
-        pl.close()
+        plt.close()
         K = np.max([[len(cor['coordinates']) for cor in cc] for cc in coors])
         cc1 = np.nan * np.zeros(np.shape(coors) + (K,))
         cc2 = np.nan * np.zeros(np.shape(coors) + (K,))
@@ -643,7 +643,7 @@ def nb_view_patches3d(Y_r, A, C, dims, image_type='mean', Yr=None,
                                                        for m in colormap(np.arange(colormap.N))])
         cmap.high = image_neurons.max()
         coors = get_contours(A, dims, thr=thr)
-        pl.close()
+        plt.close()
         cc1 = [[(l[:, 0]) for l in n['coordinates']] for n in coors]
         cc2 = [[(l[:, 1]) for l in n['coordinates']] for n in coors]
         length = np.ravel([list(map(len, cc)) for cc in cc1])
@@ -867,9 +867,9 @@ def matrixMontage(spcomps, *args, **kwargs):
     numcomps, _, _ = spcomps.shape
     rowcols = int(np.ceil(np.sqrt(numcomps)))
     for k, comp in enumerate(spcomps):
-        pl.subplot(rowcols, rowcols, k + 1)
-        pl.imshow(comp, *args, **kwargs)
-        pl.axis('off')
+        plt.subplot(rowcols, rowcols, k + 1)
+        plt.imshow(comp, *args, **kwargs)
+        plt.axis('off')
 
 
 VIDEO_TAG = """<video controls>
@@ -889,7 +889,7 @@ def anim_to_html(anim, fps=20):
     return VIDEO_TAG.format(anim._encoded_video.decode('ascii'))
 
 def display_animation(anim, fps=20):
-    pl.close(anim._fig)
+    plt.close(anim._fig)
     return HTML(anim_to_html(anim, fps=fps))
 
 def view_patches_bar(Yr, A, C, b, f, d1, d2, YrA=None, img=None,
@@ -932,7 +932,7 @@ def view_patches_bar(Yr, A, C, b, f, d1, d2, YrA=None, img=None,
             prediction values from the CNN classifier
     """
 
-    pl.ion()
+    plt.ion()
     if 'csc_matrix' not in str(type(A)):
         A = csc_matrix(A)
 
@@ -949,13 +949,13 @@ def view_patches_bar(Yr, A, C, b, f, d1, d2, YrA=None, img=None,
     if img is None:
         img = np.reshape(np.array(A.mean(axis=1)), (d1, d2), order='F')
 
-    fig = pl.figure(figsize=(10, 10))
+    fig = plt.figure(figsize=(10, 10))
 
-    axcomp = pl.axes([0.05, 0.05, 0.9, 0.03])
+    axcomp = plt.axes([0.05, 0.05, 0.9, 0.03])
 
-    ax1 = pl.axes([0.05, 0.55, 0.4, 0.4])
-    ax3 = pl.axes([0.55, 0.55, 0.4, 0.4])
-    ax2 = pl.axes([0.05, 0.1, 0.9, 0.4])
+    ax1 = plt.axes([0.05, 0.55, 0.4, 0.4])
+    ax3 = plt.axes([0.55, 0.55, 0.4, 0.4])
+    ax2 = plt.axes([0.05, 0.1, 0.9, 0.4])
 
     s_comp = matplotlib.widgets.Slider(axcomp, 'Component', 0, nr + nb - 1, valinit=0)
     vmax = np.percentile(img, 95)
@@ -968,7 +968,7 @@ def update(val):
 
             ax1.cla()
             imgtmp = np.reshape(A[:, i].toarray(), (d1, d2), order='F')
-            ax1.imshow(imgtmp, interpolation='None', cmap=pl.cm.gray, vmax=np.max(imgtmp)*0.5)
+            ax1.imshow(imgtmp, interpolation='None', cmap=matplotlib.cm.gray, vmax=np.max(imgtmp)*0.5)
             ax1.set_title('Spatial component ' + str(i + 1))
             ax1.axis('off')
 
@@ -985,11 +985,11 @@ def update(val):
                         bbox=dict(edgecolor='k', facecolor='w', alpha=.5))
 
             ax3.cla()
-            ax3.imshow(img, interpolation='None', cmap=pl.cm.gray, vmax=vmax)
+            ax3.imshow(img, interpolation='None', cmap=matplotlib.cm.gray, vmax=vmax)
             imgtmp2 = imgtmp.copy()
             imgtmp2[imgtmp2 == 0] = np.nan
             ax3.imshow(imgtmp2, interpolation='None',
-                       alpha=0.5, cmap=pl.cm.hot)
+                       alpha=0.5, cmap=matplotlib.cm.hot)
             ax3.axis('off')
         else:
             ax1.cla()
@@ -1021,7 +1021,7 @@ def arrow_key_image_control(event):
     s_comp.on_changed(update)
     s_comp.set_val(0)
     fig.canvas.mpl_connect('key_release_event', arrow_key_image_control)
-    pl.show()
+    plt.show()
 
 def plot_contours(A, Cn, thr=None, thr_method='max', maxthr=0.2, nrgthr=0.9, display_numbers=True, max_number=None,
                   cmap=None, swap_dim=False, colors='w', vmin=None, vmax=None, coordinates=None,
@@ -1081,13 +1081,13 @@ def plot_contours(A, Cn, thr=None, thr_method='max', maxthr=0.2, nrgthr=0.9, dis
             color = kwargs[key]
             kwargs.pop(key)
 
-    ax = pl.gca()
+    ax = plt.gca()
     if vmax is None and vmin is None:
-        pl.imshow(Cn, interpolation=None, cmap=cmap,
+        plt.imshow(Cn, interpolation=None, cmap=cmap,
                   vmin=np.percentile(Cn[~np.isnan(Cn)], 1),
                   vmax=np.percentile(Cn[~np.isnan(Cn)], 99))
     else:
-        pl.imshow(Cn, interpolation=None, cmap=cmap, vmin=vmin, vmax=vmax)
+        plt.imshow(Cn, interpolation=None, cmap=cmap, vmin=vmin, vmax=vmax)
 
     if coordinates is None:
         coordinates = get_contours(A, np.shape(Cn), thr, thr_method, swap_dim)
@@ -1095,7 +1095,7 @@ def plot_contours(A, Cn, thr=None, thr_method='max', maxthr=0.2, nrgthr=0.9, dis
         v = c['coordinates']
         c['bbox'] = [np.floor(np.nanmin(v[:, 1])), np.ceil(np.nanmax(v[:, 1])),
                      np.floor(np.nanmin(v[:, 0])), np.ceil(np.nanmax(v[:, 0]))]
-        pl.plot(*v.T, c=colors, **contour_args)
+        plt.plot(*v.T, c=colors, **contour_args)
 
     if display_numbers:
         d1, d2 = np.shape(Cn)
@@ -1123,15 +1123,15 @@ def GetBox(centers, R, dims):
 
     nx = int(sqrt(num_comps) * 1.3) if comps_per_row is None else comps_per_row
     ny = int(ceil(num_comps / float(nx)))
-    pl.figure(figsize=(nx, ny))
+    plt.figure(figsize=(nx, ny))
     for i, a in enumerate(Ab.T[:num_comps]):
-        ax = pl.subplot(ny, nx, i + 1)
+        ax = plt.subplot(ny, nx, i + 1)
         s = a.toarray().reshape(dims, order='F')
         box = GetBox(np.array(center_of_mass(s), dtype=np.int16), size, dims)
-        pl.imshow(smoother(s[list(map(lambda a: slice(*a), box))]),
+        plt.imshow(smoother(s[list(map(lambda a: slice(*a), box))]),
                   cmap=cmap, interpolation='nearest')
         ax.axis('off')
-    pl.subplots_adjust(0, 0, 1, 1, .06, .06)
+    plt.subplots_adjust(0, 0, 1, 1, .06, .06)
 
 
 def nb_inspect_correlation_pnr(corr, pnr, cmap='jet', num_bins=100):
@@ -1192,23 +1192,23 @@ def inspect_correlation_pnr(correlation_image_pnr, pnr_image):
             peak-to-noise image created with caiman.summary_images.correlation_pnr
     """
 
-    fig = pl.figure(figsize=(10, 4))
-    pl.axes([0.05, 0.2, 0.4, 0.7])
-    im_cn = pl.imshow(correlation_image_pnr, cmap='jet')
-    pl.title('correlation image')
-    pl.colorbar()
-    pl.axes([0.5, 0.2, 0.4, 0.7])
-    im_pnr = pl.imshow(pnr_image, cmap='jet')
-    pl.title('PNR')
-    pl.colorbar()
-
-    s_cn_max = matplotlib.widgets.Slider(pl.axes([0.05, 0.01, 0.35, 0.03]), 'vmax',
+    fig = plt.figure(figsize=(10, 4))
+    plt.axes([0.05, 0.2, 0.4, 0.7])
+    im_cn = plt.imshow(correlation_image_pnr, cmap='jet')
+    plt.title('correlation image')
+    plt.colorbar()
+    plt.axes([0.5, 0.2, 0.4, 0.7])
+    im_pnr = plt.imshow(pnr_image, cmap='jet')
+    plt.title('PNR')
+    plt.colorbar()
+
+    s_cn_max = matplotlib.widgets.Slider(plt.axes([0.05, 0.01, 0.35, 0.03]), 'vmax',
                       correlation_image_pnr.min(), correlation_image_pnr.max(), valinit=correlation_image_pnr.max())
-    s_cn_min = matplotlib.widgets.Slider(pl.axes([0.05, 0.07, 0.35, 0.03]), 'vmin',
+    s_cn_min = matplotlib.widgets.Slider(plt.axes([0.05, 0.07, 0.35, 0.03]), 'vmin',
                       correlation_image_pnr.min(), correlation_image_pnr.max(), valinit=correlation_image_pnr.min())
-    s_pnr_max = matplotlib.widgets.Slider(pl.axes([0.5, 0.01, 0.35, 0.03]), 'vmax',
+    s_pnr_max = matplotlib.widgets.Slider(plt.axes([0.5, 0.01, 0.35, 0.03]), 'vmax',
                        pnr_image.min(), pnr_image.max(), valinit=pnr_image.max())
-    s_pnr_min = matplotlib.widgets.Slider(pl.axes([0.5, 0.07, 0.35, 0.03]), 'vmin',
+    s_pnr_min = matplotlib.widgets.Slider(plt.axes([0.5, 0.07, 0.35, 0.03]), 'vmin',
                        pnr_image.min(), pnr_image.max(), valinit=pnr_image.min())
 
     def update(val):
@@ -1289,7 +1289,7 @@ def rect_draw(row_minmax: ArrayLike,
         rect: matplotlib Rectangle object
     """
     if ax is None:
-        ax = pl.gca()
+        ax = plt.gca()
         
     box_origin = (col_minmax[0], row_minmax[0])
     box_height = row_minmax[1] - row_minmax[0] 
@@ -1350,7 +1350,7 @@ def view_quilt(template_image: np.ndarray,
     patch_rows, patch_cols = get_rectangle_coords(im_dims, stride, overlap)
 
     if ax is None:             
-        f, ax = pl.subplots(figsize=figsize)
+        f, ax = plt.subplots(figsize=figsize)
         
     ax.imshow(template_image, cmap='gray', vmin=vmin, vmax=vmax)
     for patch_row in patch_rows:
diff --git a/demos/general/demo_behavior.py b/demos/general/demo_behavior.py
index 643428f78..9547d4e39 100755
--- a/demos/general/demo_behavior.py
+++ b/demos/general/demo_behavior.py
@@ -11,7 +11,7 @@
 from IPython import get_ipython
 import logging
 import numpy as np
-import pylab as pl
+import matplotlib.pyplot as plt
 
 try:
     cv2.setNumThreads(0)
@@ -47,7 +47,7 @@ def main():
     pass # For compatibility between running under IDE and the CLI
 
     #%%
-    pl.ion()
+    plt.ion()
 
     fname = [u'demo_behavior.h5']
     if fname[0] in ['demo_behavior.h5']:
@@ -87,9 +87,9 @@ def main():
     
     #%%
     idd = 0
-    axlin = pl.subplot(n_components, 2, 2)
+    axlin = plt.subplot(n_components, 2, 2)
     for mag, dirct, spatial_filter in zip(mags, dircts_thresh, spatial_filter_):
-        pl.subplot(n_components, 2, 1 + idd * 2)
+        plt.subplot(n_components, 2, 1 + idd * 2)
         min_x, min_y = np.min(np.where(mask), 1)
 
         spfl = spatial_filter
@@ -99,14 +99,14 @@ def main():
 
         mask[min_x:max_x, min_y:max_y] = spfl
         mask[mask < np.nanpercentile(spfl, 70)] = np.nan
-        pl.imshow(m[0], cmap='gray')
-        pl.imshow(mask, alpha=.5)
-        pl.axis('off')
+        plt.imshow(m[0], cmap='gray')
+        plt.imshow(mask, alpha=.5)
+        plt.axis('off')
 
-        axelin = pl.subplot(n_components, 2, 2 + idd * 2, sharex=axlin)
-        pl.plot(mag / 10, 'k')
+        axelin = plt.subplot(n_components, 2, 2 + idd * 2, sharex=axlin)
+        plt.plot(mag / 10, 'k')
         dirct[mag < 0.5 * np.std(mag)] = np.nan
-        pl.plot(dirct, 'r-', linewidth=2)
+        plt.plot(dirct, 'r-', linewidth=2)
 
         idd += 1