diff --git a/deeplabcut/create_project/add.py b/deeplabcut/create_project/add.py
index c79a216e2e..3fa68ab5dc 100644
--- a/deeplabcut/create_project/add.py
+++ b/deeplabcut/create_project/add.py
@@ -30,7 +30,7 @@ def add_new_videos(config, videos, copy_videos=False, coords=None):
     --------
     Video will be added, with cropping dimenions according to the frame dimensinos of mouse5.avi
     >>> deeplabcut.add_new_videos('/home/project/reaching-task-Tanmay-2018-08-23/config.yaml',['/data/videos/mouse5.avi'])
-    
+
     Video will be added, with cropping dimenions [0,100,0,200]
     >>> deeplabcut.add_new_videos('/home/project/reaching-task-Tanmay-2018-08-23/config.yaml',['/data/videos/mouse5.avi'],copy_videos=False,coords=[[0,100,0,200]])
 
diff --git a/deeplabcut/create_project/demo_data.py b/deeplabcut/create_project/demo_data.py
index 32e5c31a8d..176c8ba7fc 100644
--- a/deeplabcut/create_project/demo_data.py
+++ b/deeplabcut/create_project/demo_data.py
@@ -17,18 +17,22 @@
 
 def load_demo_data(config, createtrainingset=True):
     """
-  Loads the demo data. Make sure that you are in the same directory where you have downloaded or cloned the deeplabcut.
-
-  Parameter
-    ----------
-    config : string
-        Full path of the config.yaml file of the provided demo dataset as a string.
-
-    Example
-    --------
-    >>> deeplabcut.load_demo_data('config.yaml')
-    --------
-  """
+    Loads the demo data -- subset from trail-tracking data in Mathis et al. 2018. 
+    When loading, it sets paths correctly to run this project on your system
+    
+    Parameter
+      ----------
+      config : string
+          Full path of the config.yaml file of the provided demo dataset as a string.
+
+      createtrainingset : bool 
+          Boolean variable indicating if a training set shall be created. 
+          
+      Example
+      --------
+      >>> deeplabcut.load_demo_data('config.yaml')
+      --------
+    """
     config = Path(config).resolve()
     config = str(config)
 
diff --git a/deeplabcut/generate_training_dataset/multiple_individuals_trainingsetmanipulation.py b/deeplabcut/generate_training_dataset/multiple_individuals_trainingsetmanipulation.py
index 69442301e6..822fc0ea7a 100755
--- a/deeplabcut/generate_training_dataset/multiple_individuals_trainingsetmanipulation.py
+++ b/deeplabcut/generate_training_dataset/multiple_individuals_trainingsetmanipulation.py
@@ -31,7 +31,10 @@
 
 
 def format_multianimal_training_data(
-    df, train_inds, project_path, n_decimals=2,
+    df,
+    train_inds,
+    project_path,
+    n_decimals=2,
 ):
     train_data = []
     nrows = df.shape[0]
@@ -283,7 +286,10 @@ def create_multianimaltraining_dataset(
 
         # Make training file!
         data = format_multianimal_training_data(
-            Data, trainIndices, cfg["project_path"], numdigits,
+            Data,
+            trainIndices,
+            cfg["project_path"],
+            numdigits,
         )
 
         if len(trainIndices) > 0:
@@ -339,7 +345,10 @@ def create_multianimaltraining_dataset(
             )
             path_test_config = str(
                 os.path.join(
-                    cfg["project_path"], Path(modelfoldername), "test", "pose_cfg.yaml",
+                    cfg["project_path"],
+                    Path(modelfoldername),
+                    "test",
+                    "pose_cfg.yaml",
                 )
             )
             path_inference_config = str(
@@ -440,7 +449,10 @@ def create_multianimaltraining_dataset(
 
 
 def convert_cropped_to_standard_dataset(
-    config_path, recreate_datasets=True, delete_crops=True, back_up=True,
+    config_path,
+    recreate_datasets=True,
+    delete_crops=True,
+    back_up=True,
 ):
     import pandas as pd
     import pickle
@@ -479,7 +491,8 @@ def convert_cropped_to_standard_dataset(
         return
 
     datasets_folder = os.path.join(
-        project_path, auxiliaryfunctions.GetTrainingSetFolder(cfg),
+        project_path,
+        auxiliaryfunctions.GetTrainingSetFolder(cfg),
     )
     df_old = pd.read_hdf(
         os.path.join(datasets_folder, "CollectedData_" + cfg["scorer"] + ".h5"),
diff --git a/deeplabcut/generate_training_dataset/trainingsetmanipulation.py b/deeplabcut/generate_training_dataset/trainingsetmanipulation.py
index 0b43dad840..b02d365250 100755
--- a/deeplabcut/generate_training_dataset/trainingsetmanipulation.py
+++ b/deeplabcut/generate_training_dataset/trainingsetmanipulation.py
@@ -318,7 +318,7 @@ def check_labels(
 
 
 def boxitintoacell(joints):
-    """ Auxiliary function for creating matfile."""
+    """Auxiliary function for creating matfile."""
     outer = np.array([[None]], dtype=object)
     outer[0, 0] = np.array(joints, dtype="int64")
     return outer
@@ -460,9 +460,11 @@ def merge_annotateddatasets(cfg, trainingsetfolder_full):
 
 
 def SplitTrials(
-    trialindex, trainFraction=0.8, enforce_train_fraction=False,
+    trialindex,
+    trainFraction=0.8,
+    enforce_train_fraction=False,
 ):
-    """ Split a trial index into train and test sets. Also checks that the trainFraction is a two digit number between 0 an 1. The reason
+    """Split a trial index into train and test sets. Also checks that the trainFraction is a two digit number between 0 an 1. The reason
     is that the folders contain the trainfraction as int(100*trainFraction).
     If enforce_train_fraction is True, train and test indices are padded with -1
     such that the ratio of their lengths is exactly the desired train fraction.
@@ -487,7 +489,9 @@ def SplitTrials(
         train_indices = shuffle[: int(train_size)]
         if enforce_train_fraction and not train_size.is_integer():
             train_indices, test_indices = pad_train_test_indices(
-                train_indices, test_indices, train_fraction,
+                train_indices,
+                test_indices,
+                train_fraction,
             )
         return train_indices, test_indices
 
@@ -505,7 +509,8 @@ def pad_train_test_indices(train_inds, test_inds, train_fraction):
     min_n_train = int(round(min_length_req * train_fraction))
     min_n_test = min_length_req - min_n_train
     mult = max(
-        math.ceil(n_train_inds / min_n_train), math.ceil(n_test_inds / min_n_test),
+        math.ceil(n_train_inds / min_n_train),
+        math.ceil(n_test_inds / min_n_test),
     )
     n_train = mult * min_n_train
     n_test = mult * min_n_test
@@ -570,7 +575,8 @@ def mergeandsplit(config, trainindex=0, uniform=True):
         Data = pd.read_hdf(fn + ".h5")
     except FileNotFoundError:
         Data = merge_annotateddatasets(
-            cfg, Path(os.path.join(project_path, trainingsetfolder)),
+            cfg,
+            Path(os.path.join(project_path, trainingsetfolder)),
         )
         if Data is None:
             return [], []
@@ -582,7 +588,9 @@ def mergeandsplit(config, trainindex=0, uniform=True):
         TrainingFraction = cfg["TrainingFraction"]
         trainFraction = TrainingFraction[trainindex]
         trainIndices, testIndices = SplitTrials(
-            range(len(Data.index)), trainFraction, True,
+            range(len(Data.index)),
+            trainFraction,
+            True,
         )
     else:  # leave one folder out split
         videos = cfg["video_sets"].keys()
@@ -733,7 +741,8 @@ def create_training_dataset(
         )
 
         Data = merge_annotateddatasets(
-            cfg, Path(os.path.join(project_path, trainingsetfolder)),
+            cfg,
+            Path(os.path.join(project_path, trainingsetfolder)),
         )
         if Data is None:
             return
@@ -946,19 +955,21 @@ def create_training_dataset(
 
 
 def get_largestshuffle_index(config):
-    """ Returns the largest shuffle for all dlc-models in the current iteration."""
+    """Returns the largest shuffle for all dlc-models in the current iteration."""
     cfg = auxiliaryfunctions.read_config(config)
     project_path = cfg["project_path"]
     iterate = "iteration-" + str(cfg["iteration"])
     dlc_model_path = os.path.join(project_path, "dlc-models", iterate)
     if os.path.isdir(dlc_model_path):
         models = os.listdir(dlc_model_path)
-        # sort the models directories
+        # sort the model directories
         models.sort(key=lambda f: int("".join(filter(str.isdigit, f))))
-        # get the shuffle index
-        max_shuffle_index = int(models[-1].split("shuffle")[-1])
+
+        # get the shuffle index and offset by 1.
+        max_shuffle_index = int(models[-1].split("shuffle")[-1]) + 1
     else:
         max_shuffle_index = 0
+
     return max_shuffle_index
 
 
@@ -967,7 +978,7 @@ def create_training_model_comparison(
     trainindex=0,
     num_shuffles=1,
     net_types=["resnet_50"],
-    augmenter_types=["default"],
+    augmenter_types=["imgaug"],
     userfeedback=False,
     windows2linux=False,
 ):
@@ -1001,12 +1012,19 @@ def create_training_model_comparison(
         If this is set to false, then all requested train/test splits are created (no matter if they already exist). If you
         want to assure that previous splits etc. are not overwritten, then set this to True and you will be asked for each split.
 
+    Returns
+    ----------
+    shuffle_list: list
+        List of indices corresponding to the trainigsplits/models that were created.
+
     Example
     --------
-    >>> deeplabcut.create_training_model_comparison('/analysis/project/reaching-task/config.yaml',num_shuffles=1,net_types=['resnet_50','resnet_152'],augmenter_types=['tensorpack','deterministic'])
+    >>> shuffle_list = deeplabcut.create_training_model_comparison('/analysis/project/reaching-task/config.yaml',num_shuffles=1,net_types=['resnet_50','resnet_152'],augmenter_types=['tensorpack','deterministic'])
 
     Windows:
-    >>> deeplabcut.create_training_model_comparison('C:\\Users\\Ulf\\looming-task\\config.yaml',num_shuffles=1,net_types=['resnet_50','resnet_152'],augmenter_types=['tensorpack','deterministic'])
+    >>> shuffle_list = deeplabcut.create_training_model_comparison('C:\\Users\\Ulf\\looming-task\\config.yaml',num_shuffles=1,net_types=['resnet_50','resnet_152'],augmenter_types=['tensorpack','deterministic'])
+
+    See examples/testscript_openfielddata_augmentationcomparison.py for an example of how to use shuffle_list.
 
     --------
     """
@@ -1034,6 +1052,7 @@ def create_training_model_comparison(
 
     largestshuffleindex = get_largestshuffle_index(config)
 
+    shuffle_list = []
     for shuffle in range(num_shuffles):
         trainIndices, testIndices = mergeandsplit(
             config, trainindex=trainindex, uniform=True
@@ -1046,6 +1065,8 @@ def create_training_model_comparison(
                     + idx_net * len(augmenter_types)
                     + shuffle * len(augmenter_types) * len(net_types)
                 )
+
+                shuffle_list.append(get_max_shuffle_idx)
                 log_info = str(
                     "Shuffle index:"
                     + str(get_max_shuffle_idx)
@@ -1055,6 +1076,8 @@ def create_training_model_comparison(
                     + aug
                     + ", trainsetindex:"
                     + str(trainindex)
+                    + ", frozen shuffle ID:"
+                    + str(shuffle)
                 )
                 create_training_dataset(
                     config,
@@ -1066,3 +1089,5 @@ def create_training_model_comparison(
                     userfeedback=userfeedback,
                 )
                 logger.info(log_info)
+
+    return shuffle_list
diff --git a/deeplabcut/gui/outlier_frame_extraction_toolbox.py b/deeplabcut/gui/outlier_frame_extraction_toolbox.py
index 393bd8a031..e48bbe8491 100644
--- a/deeplabcut/gui/outlier_frame_extraction_toolbox.py
+++ b/deeplabcut/gui/outlier_frame_extraction_toolbox.py
@@ -42,7 +42,6 @@ def getColorIndices(self, img, bodyparts):
         """
         Returns the colormaps ticks and . The order of ticks labels is reversed.
         """
-        #        im = io.imread(img)
         norm = mcolors.Normalize(vmin=np.min(img), vmax=np.max(img))
         ticks = np.linspace(np.min(img), np.max(img), len(bodyparts))[::-1]
         return norm, ticks
diff --git a/deeplabcut/pose_estimation_tensorflow/backbones/efficientnet_model.py b/deeplabcut/pose_estimation_tensorflow/backbones/efficientnet_model.py
index 1f82dd62d1..ba7a41c23a 100644
--- a/deeplabcut/pose_estimation_tensorflow/backbones/efficientnet_model.py
+++ b/deeplabcut/pose_estimation_tensorflow/backbones/efficientnet_model.py
@@ -386,7 +386,7 @@ def call(
 
 class Model(tf.keras.Model):
     """A class implements tf.keras.Model for MNAS-like model.
-      Reference: https://arxiv.org/abs/1807.11626
+    Reference: https://arxiv.org/abs/1807.11626
     """
 
     def __init__(self, blocks_args=None, global_params=None):
diff --git a/deeplabcut/pose_estimation_tensorflow/backbones/mobilenet.py b/deeplabcut/pose_estimation_tensorflow/backbones/mobilenet.py
index 2b642bcc81..a1c911e1e9 100644
--- a/deeplabcut/pose_estimation_tensorflow/backbones/mobilenet.py
+++ b/deeplabcut/pose_estimation_tensorflow/backbones/mobilenet.py
@@ -35,13 +35,13 @@ def apply_activation(x, name=None, activation_fn=None):
 def _set_arg_scope_defaults(defaults):
     """Sets arg scope defaults for all items present in defaults.
 
-  Args:
-    defaults: dictionary/list of pairs, containing a mapping from
-    function to a dictionary of default args.
+    Args:
+      defaults: dictionary/list of pairs, containing a mapping from
+      function to a dictionary of default args.
 
-  Yields:
-    context manager where all defaults are set.
-  """
+    Yields:
+      context manager where all defaults are set.
+    """
     if hasattr(defaults, "items"):
         items = list(defaults.items())
     else:
@@ -88,16 +88,16 @@ def __exit__(self, exc_type, exc_value, traceback):
 def safe_arg_scope(funcs, **kwargs):
     """Returns `slim.arg_scope` with all None arguments removed.
 
-  Arguments:
-    funcs: Functions to pass to `arg_scope`.
-    **kwargs: Arguments to pass to `arg_scope`.
+    Arguments:
+      funcs: Functions to pass to `arg_scope`.
+      **kwargs: Arguments to pass to `arg_scope`.
 
-  Returns:
-    arg_scope or No-op context manager.
+    Returns:
+      arg_scope or No-op context manager.
 
-  Note: can be useful if None value should be interpreted as "do not overwrite
-    this parameter value".
-  """
+    Note: can be useful if None value should be interpreted as "do not overwrite
+      this parameter value".
+    """
     filtered_args = {name: value for name, value in kwargs.items() if value is not None}
     if filtered_args:
         return slim.arg_scope(funcs, **filtered_args)
@@ -118,55 +118,55 @@ def mobilenet_base(  # pylint: disable=invalid-name
 ):
     """Mobilenet base network.
 
-  Constructs a network from inputs to the given final endpoint. By default
-  the network is constructed in inference mode. To create network
-  in training mode use:
-
-  with slim.arg_scope(mobilenet.training_scope()):
-     logits, endpoints = mobilenet_base(...)
-
-  Args:
-    inputs: a tensor of shape [batch_size, height, width, channels].
-    conv_defs: A list of op(...) layers specifying the net architecture.
-    multiplier: Float multiplier for the depth (number of channels)
-      for all convolution ops. The value must be greater than zero. Typical
-      usage will be to set this value in (0, 1) to reduce the number of
-      parameters or computation cost of the model.
-    final_endpoint: The name of last layer, for early termination for
-    for V1-based networks: last layer is "layer_14", for V2: "layer_20"
-    output_stride: An integer that specifies the requested ratio of input to
-      output spatial resolution. If not None, then we invoke atrous convolution
-      if necessary to prevent the network from reducing the spatial resolution
-      of the activation maps. Allowed values are 1 or any even number, excluding
-      zero. Typical values are 8 (accurate fully convolutional mode), 16
-      (fast fully convolutional mode), and 32 (classification mode).
-
-      NOTE- output_stride relies on all consequent operators to support dilated
-      operators via "rate" parameter. This might require wrapping non-conv
-      operators to operate properly.
-
-    use_explicit_padding: Use 'VALID' padding for convolutions, but prepad
-      inputs so that the output dimensions are the same as if 'SAME' padding
-      were used.
-    scope: optional variable scope.
-    is_training: How to setup batch_norm and other ops. Note: most of the time
-      this does not need be set directly. Use mobilenet.training_scope() to set
-      up training instead. This parameter is here for backward compatibility
-      only. It is safe to set it to the value matching
-      training_scope(is_training=...). It is also safe to explicitly set
-      it to False, even if there is outer training_scope set to to training.
-      (The network will be built in inference mode). If this is set to None,
-      no arg_scope is added for slim.batch_norm's is_training parameter.
-
-  Returns:
-    tensor_out: output tensor.
-    end_points: a set of activations for external use, for example summaries or
-                losses.
-
-  Raises:
-    ValueError: depth_multiplier <= 0, or the target output_stride is not
-                allowed.
-  """
+    Constructs a network from inputs to the given final endpoint. By default
+    the network is constructed in inference mode. To create network
+    in training mode use:
+
+    with slim.arg_scope(mobilenet.training_scope()):
+       logits, endpoints = mobilenet_base(...)
+
+    Args:
+      inputs: a tensor of shape [batch_size, height, width, channels].
+      conv_defs: A list of op(...) layers specifying the net architecture.
+      multiplier: Float multiplier for the depth (number of channels)
+        for all convolution ops. The value must be greater than zero. Typical
+        usage will be to set this value in (0, 1) to reduce the number of
+        parameters or computation cost of the model.
+      final_endpoint: The name of last layer, for early termination for
+      for V1-based networks: last layer is "layer_14", for V2: "layer_20"
+      output_stride: An integer that specifies the requested ratio of input to
+        output spatial resolution. If not None, then we invoke atrous convolution
+        if necessary to prevent the network from reducing the spatial resolution
+        of the activation maps. Allowed values are 1 or any even number, excluding
+        zero. Typical values are 8 (accurate fully convolutional mode), 16
+        (fast fully convolutional mode), and 32 (classification mode).
+
+        NOTE- output_stride relies on all consequent operators to support dilated
+        operators via "rate" parameter. This might require wrapping non-conv
+        operators to operate properly.
+
+      use_explicit_padding: Use 'VALID' padding for convolutions, but prepad
+        inputs so that the output dimensions are the same as if 'SAME' padding
+        were used.
+      scope: optional variable scope.
+      is_training: How to setup batch_norm and other ops. Note: most of the time
+        this does not need be set directly. Use mobilenet.training_scope() to set
+        up training instead. This parameter is here for backward compatibility
+        only. It is safe to set it to the value matching
+        training_scope(is_training=...). It is also safe to explicitly set
+        it to False, even if there is outer training_scope set to to training.
+        (The network will be built in inference mode). If this is set to None,
+        no arg_scope is added for slim.batch_norm's is_training parameter.
+
+    Returns:
+      tensor_out: output tensor.
+      end_points: a set of activations for external use, for example summaries or
+                  losses.
+
+    Raises:
+      ValueError: depth_multiplier <= 0, or the target output_stride is not
+                  allowed.
+    """
     if multiplier <= 0:
         raise ValueError("multiplier is not greater than zero.")
 
@@ -283,43 +283,43 @@ def mobilenet(
 ):
     """Mobilenet model for classification, supports both V1 and V2.
 
-  Note: default mode is inference, use mobilenet.training_scope to create
-  training network.
-
-
-  Args:
-    inputs: a tensor of shape [batch_size, height, width, channels].
-    num_classes: number of predicted classes. If 0 or None, the logits layer
-      is omitted and the input features to the logits layer (before dropout)
-      are returned instead.
-    prediction_fn: a function to get predictions out of logits
-      (default softmax).
-    reuse: whether or not the network and its variables should be reused. To be
-      able to reuse 'scope' must be given.
-    scope: Optional variable_scope.
-    base_only: if True will only create the base of the network (no pooling
-    and no logits).
-    **mobilenet_args: passed to mobilenet_base verbatim.
-      - conv_defs: list of conv defs
-      - multiplier: Float multiplier for the depth (number of channels)
-      for all convolution ops. The value must be greater than zero. Typical
-      usage will be to set this value in (0, 1) to reduce the number of
-      parameters or computation cost of the model.
-      - output_stride: will ensure that the last layer has at most total stride.
-      If the architecture calls for more stride than that provided
-      (e.g. output_stride=16, but the architecture has 5 stride=2 operators),
-      it will replace output_stride with fractional convolutions using Atrous
-      Convolutions.
-
-  Returns:
-    logits: the pre-softmax activations, a tensor of size
-      [batch_size, num_classes]
-    end_points: a dictionary from components of the network to the corresponding
-      activation tensor.
-
-  Raises:
-    ValueError: Input rank is invalid.
-  """
+    Note: default mode is inference, use mobilenet.training_scope to create
+    training network.
+
+
+    Args:
+      inputs: a tensor of shape [batch_size, height, width, channels].
+      num_classes: number of predicted classes. If 0 or None, the logits layer
+        is omitted and the input features to the logits layer (before dropout)
+        are returned instead.
+      prediction_fn: a function to get predictions out of logits
+        (default softmax).
+      reuse: whether or not the network and its variables should be reused. To be
+        able to reuse 'scope' must be given.
+      scope: Optional variable_scope.
+      base_only: if True will only create the base of the network (no pooling
+      and no logits).
+      **mobilenet_args: passed to mobilenet_base verbatim.
+        - conv_defs: list of conv defs
+        - multiplier: Float multiplier for the depth (number of channels)
+        for all convolution ops. The value must be greater than zero. Typical
+        usage will be to set this value in (0, 1) to reduce the number of
+        parameters or computation cost of the model.
+        - output_stride: will ensure that the last layer has at most total stride.
+        If the architecture calls for more stride than that provided
+        (e.g. output_stride=16, but the architecture has 5 stride=2 operators),
+        it will replace output_stride with fractional convolutions using Atrous
+        Convolutions.
+
+    Returns:
+      logits: the pre-softmax activations, a tensor of size
+        [batch_size, num_classes]
+      end_points: a dictionary from components of the network to the corresponding
+        activation tensor.
+
+    Raises:
+      ValueError: Input rank is invalid.
+    """
     is_training = mobilenet_args.get("is_training", False)
     input_shape = inputs.get_shape().as_list()
     if len(input_shape) != 4:
@@ -363,15 +363,15 @@ def mobilenet(
 def global_pool(input_tensor, pool_op=tf.nn.avg_pool2d):
     """Applies avg pool to produce 1x1 output.
 
-  NOTE: This function is funcitonally equivalenet to reduce_mean, but it has
-  baked in average pool which has better support across hardware.
+    NOTE: This function is funcitonally equivalenet to reduce_mean, but it has
+    baked in average pool which has better support across hardware.
 
-  Args:
-    input_tensor: input tensor
-    pool_op: pooling op (avg pool is default)
-  Returns:
-    a tensor batch_size x 1 x 1 x depth.
-  """
+    Args:
+      input_tensor: input tensor
+      pool_op: pooling op (avg pool is default)
+    Returns:
+      a tensor batch_size x 1 x 1 x depth.
+    """
     shape = input_tensor.get_shape().as_list()
     if shape[1] is None or shape[2] is None:
         kernel_size = tf.convert_to_tensor(
@@ -401,28 +401,28 @@ def training_scope(
 ):
     """Defines Mobilenet training scope.
 
-  Usage:
-     with tf.contrib.slim.arg_scope(mobilenet.training_scope()):
-       logits, endpoints = mobilenet_v2.mobilenet(input_tensor)
-
-     # the network created will be trainble with dropout/batch norm
-     # initialized appropriately.
-  Args:
-    is_training: if set to False this will ensure that all customizations are
-      set to non-training mode. This might be helpful for code that is reused
-      across both training/evaluation, but most of the time training_scope with
-      value False is not needed. If this is set to None, the parameters is not
-      added to the batch_norm arg_scope.
-
-    weight_decay: The weight decay to use for regularizing the model.
-    stddev: Standard deviation for initialization, if negative uses xavier.
-    dropout_keep_prob: dropout keep probability (not set if equals to None).
-    bn_decay: decay for the batch norm moving averages (not set if equals to
-      None).
-
-  Returns:
-    An argument scope to use via arg_scope.
-  """
+    Usage:
+       with tf.contrib.slim.arg_scope(mobilenet.training_scope()):
+         logits, endpoints = mobilenet_v2.mobilenet(input_tensor)
+
+       # the network created will be trainble with dropout/batch norm
+       # initialized appropriately.
+    Args:
+      is_training: if set to False this will ensure that all customizations are
+        set to non-training mode. This might be helpful for code that is reused
+        across both training/evaluation, but most of the time training_scope with
+        value False is not needed. If this is set to None, the parameters is not
+        added to the batch_norm arg_scope.
+
+      weight_decay: The weight decay to use for regularizing the model.
+      stddev: Standard deviation for initialization, if negative uses xavier.
+      dropout_keep_prob: dropout keep probability (not set if equals to None).
+      bn_decay: decay for the batch norm moving averages (not set if equals to
+        None).
+
+    Returns:
+      An argument scope to use via arg_scope.
+    """
     # Note: do not introduce parameters that would change the inference
     # model here (for example whether to use bias), modify conv_def instead.
     batch_norm_params = {"decay": bn_decay, "is_training": is_training}
diff --git a/deeplabcut/pose_estimation_tensorflow/backbones/mobilenet_v2.py b/deeplabcut/pose_estimation_tensorflow/backbones/mobilenet_v2.py
index 4d2aed38d7..c1cc3d7b09 100644
--- a/deeplabcut/pose_estimation_tensorflow/backbones/mobilenet_v2.py
+++ b/deeplabcut/pose_estimation_tensorflow/backbones/mobilenet_v2.py
@@ -101,40 +101,40 @@ def mobilenet(
 ):
     """Creates mobilenet V2 network.
 
-  Inference mode is created by default. To create training use training_scope
-  below.
-
-  with tf.contrib.slim.arg_scope(mobilenet_v2.training_scope()):
-     logits, endpoints = mobilenet_v2.mobilenet(input_tensor)
-
-  Args:
-    input_tensor: The input tensor
-    num_classes: number of classes
-    depth_multiplier: The multiplier applied to scale number of
-    channels in each layer.
-    scope: Scope of the operator
-    conv_defs: Allows to override default conv def.
-    finegrain_classification_mode: When set to True, the model
-    will keep the last layer large even for small multipliers. Following
-    https://arxiv.org/abs/1801.04381
-    suggests that it improves performance for ImageNet-type of problems.
-      *Note* ignored if final_endpoint makes the builder exit earlier.
-    min_depth: If provided, will ensure that all layers will have that
-    many channels after application of depth multiplier.
-    divisible_by: If provided will ensure that all layers # channels
-    will be divisible by this number.
-    activation_fn: Activation function to use, defaults to tf.nn.relu6 if not
-      specified.
-    **kwargs: passed directly to mobilenet.mobilenet:
-      prediction_fn- what prediction function to use.
-      reuse-: whether to reuse variables (if reuse set to true, scope
-      must be given).
-  Returns:
-    logits/endpoints pair
-
-  Raises:
-    ValueError: On invalid arguments
-  """
+    Inference mode is created by default. To create training use training_scope
+    below.
+
+    with tf.contrib.slim.arg_scope(mobilenet_v2.training_scope()):
+       logits, endpoints = mobilenet_v2.mobilenet(input_tensor)
+
+    Args:
+      input_tensor: The input tensor
+      num_classes: number of classes
+      depth_multiplier: The multiplier applied to scale number of
+      channels in each layer.
+      scope: Scope of the operator
+      conv_defs: Allows to override default conv def.
+      finegrain_classification_mode: When set to True, the model
+      will keep the last layer large even for small multipliers. Following
+      https://arxiv.org/abs/1801.04381
+      suggests that it improves performance for ImageNet-type of problems.
+        *Note* ignored if final_endpoint makes the builder exit earlier.
+      min_depth: If provided, will ensure that all layers will have that
+      many channels after application of depth multiplier.
+      divisible_by: If provided will ensure that all layers # channels
+      will be divisible by this number.
+      activation_fn: Activation function to use, defaults to tf.nn.relu6 if not
+        specified.
+      **kwargs: passed directly to mobilenet.mobilenet:
+        prediction_fn- what prediction function to use.
+        reuse-: whether to reuse variables (if reuse set to true, scope
+        must be given).
+    Returns:
+      logits/endpoints pair
+
+    Raises:
+      ValueError: On invalid arguments
+    """
     if conv_defs is None:
         conv_defs = V2_DEF
     if "multiplier" in kwargs:
@@ -205,23 +205,23 @@ def mobilenet_base(input_tensor, depth_multiplier=1.0, **kwargs):
 def training_scope(**kwargs):
     """Defines MobilenetV2 training scope.
 
-  Usage:
-     with tf.contrib.slim.arg_scope(mobilenet_v2.training_scope()):
-       logits, endpoints = mobilenet_v2.mobilenet(input_tensor)
+    Usage:
+       with tf.contrib.slim.arg_scope(mobilenet_v2.training_scope()):
+         logits, endpoints = mobilenet_v2.mobilenet(input_tensor)
 
-  with slim.
+    with slim.
 
-  Args:
-    **kwargs: Passed to mobilenet.training_scope. The following parameters
-    are supported:
-      weight_decay- The weight decay to use for regularizing the model.
-      stddev-  Standard deviation for initialization, if negative uses xavier.
-      dropout_keep_prob- dropout keep probability
-      bn_decay- decay for the batch norm moving averages.
+    Args:
+      **kwargs: Passed to mobilenet.training_scope. The following parameters
+      are supported:
+        weight_decay- The weight decay to use for regularizing the model.
+        stddev-  Standard deviation for initialization, if negative uses xavier.
+        dropout_keep_prob- dropout keep probability
+        bn_decay- decay for the batch norm moving averages.
 
-  Returns:
-    An `arg_scope` to use for the mobilenet v2 model.
-  """
+    Returns:
+      An `arg_scope` to use for the mobilenet v2 model.
+    """
     return lib.training_scope(**kwargs)
 
 
diff --git a/deeplabcut/pose_estimation_tensorflow/core/evaluate.py b/deeplabcut/pose_estimation_tensorflow/core/evaluate.py
index ecc61d8e5c..3c3b158144 100644
--- a/deeplabcut/pose_estimation_tensorflow/core/evaluate.py
+++ b/deeplabcut/pose_estimation_tensorflow/core/evaluate.py
@@ -18,7 +18,7 @@
 
 
 def pairwisedistances(DataCombined, scorer1, scorer2, pcutoff=-1, bodyparts=None):
-    """ Calculates the pairwise Euclidean distance metric over body parts vs. images"""
+    """Calculates the pairwise Euclidean distance metric over body parts vs. images"""
     mask = DataCombined[scorer2].xs("likelihood", level=1, axis=1) >= pcutoff
     if bodyparts == None:
         Pointwisesquareddistance = (DataCombined[scorer1] - DataCombined[scorer2]) ** 2
@@ -191,7 +191,7 @@ def calculatepafdistancebounds(
 def Plotting(
     cfg, comparisonbodyparts, DLCscorer, trainIndices, DataCombined, foldername
 ):
-    """ Function used for plotting GT and predictions """
+    """Function used for plotting GT and predictions"""
     from deeplabcut.utils import visualization
 
     colors = visualization.get_cmap(len(comparisonbodyparts), name=cfg["colormap"])
@@ -631,8 +631,10 @@ def evaluate_network(
         )
 
         # Get list of body parts to evaluate network for
-        comparisonbodyparts = auxiliaryfunctions.IntersectionofBodyPartsandOnesGivenbyUser(
-            cfg, comparisonbodyparts
+        comparisonbodyparts = (
+            auxiliaryfunctions.IntersectionofBodyPartsandOnesGivenbyUser(
+                cfg, comparisonbodyparts
+            )
         )
         # Make folder for evaluation
         auxiliaryfunctions.attempttomakefolder(
@@ -788,7 +790,7 @@ def evaluate_network(
                         for imageindex, imagename in tqdm(enumerate(Data.index)):
                             image = imread(
                                 os.path.join(cfg["project_path"], *imagename),
-                                mode="RGB",
+                                mode="skimage",
                             )
                             if scale != 1:
                                 image = imresize(image, scale)
diff --git a/deeplabcut/pose_estimation_tensorflow/core/evaluate_multianimal.py b/deeplabcut/pose_estimation_tensorflow/core/evaluate_multianimal.py
index 27374731d5..045ce1d10b 100644
--- a/deeplabcut/pose_estimation_tensorflow/core/evaluate_multianimal.py
+++ b/deeplabcut/pose_estimation_tensorflow/core/evaluate_multianimal.py
@@ -12,13 +12,9 @@
 import os
 import pickle
 from pathlib import Path
-
 import numpy as np
 import pandas as pd
-import skimage.color
 from scipy.spatial import cKDTree
-from skimage import io
-from skimage.util import img_as_ubyte
 from tqdm import tqdm
 
 from deeplabcut.pose_estimation_tensorflow.core.evaluate import make_results_file
@@ -117,6 +113,7 @@ def evaluate_multianimal_full(
     from deeplabcut.utils import (
         auxiliaryfunctions,
         auxfun_multianimal,
+        auxfun_videos,
         conversioncode,
     )
 
@@ -313,12 +310,7 @@ def evaluate_multianimal_full(
                         print("Network Evaluation underway...")
                         for imageindex, imagename in tqdm(enumerate(Data.index)):
                             image_path = os.path.join(cfg["project_path"], *imagename)
-                            image = io.imread(image_path)
-                            if image.ndim == 2 or image.shape[-1] == 1:
-                                image = skimage.color.gray2rgb(image)
-                            elif image.shape[-1] == 4:
-                                image = skimage.color.rgba2rgb(image)
-                            frame = img_as_ubyte(image)
+                            frame = auxfun_videos.imread(image_path, mode="skimage")
 
                             GT = Data.iloc[imageindex]
                             if not GT.any():
@@ -343,7 +335,8 @@ def evaluate_multianimal_full(
                             # is (sample_index, peak_y, peak_x, bpt_index) to slice the PAFs.
                             temp = df.reset_index(level="bodyparts").dropna()
                             temp["bodyparts"].replace(
-                                dict(zip(joints, range(len(joints)))), inplace=True,
+                                dict(zip(joints, range(len(joints)))),
+                                inplace=True,
                             )
                             temp["sample"] = 0
                             peaks_gt = temp.loc[
@@ -581,10 +574,8 @@ def evaluate_multianimal_full(
                         for k, v in tqdm(assemblies.items()):
                             imname = image_paths[k]
                             image_path = os.path.join(cfg["project_path"], *imname)
-                            image = io.imread(image_path)
-                            if image.ndim == 2 or image.shape[-1] == 1:
-                                image = skimage.color.gray2rgb(image)
-                            frame = img_as_ubyte(image)
+                            frame = auxfun_videos.imread(image_path, mode="skimage")
+
                             h, w, _ = np.shape(frame)
                             fig.set_size_inches(w / 100, h / 100)
                             ax.set_xlim(0, w)
@@ -619,7 +610,10 @@ def evaluate_multianimal_full(
                                 ax=ax,
                             )
                             visualization.save_labeled_frame(
-                                fig, image_path, foldername, k in trainIndices,
+                                fig,
+                                image_path,
+                                foldername,
+                                k in trainIndices,
                             )
                             visualization.erase_artists(ax)
 
diff --git a/deeplabcut/pose_estimation_tensorflow/core/predict.py b/deeplabcut/pose_estimation_tensorflow/core/predict.py
index 182cb08ac8..110ead75c2 100644
--- a/deeplabcut/pose_estimation_tensorflow/core/predict.py
+++ b/deeplabcut/pose_estimation_tensorflow/core/predict.py
@@ -61,7 +61,7 @@ def setup_pose_prediction(cfg, allow_growth=False):
 
 
 def extract_cnn_output(outputs_np, cfg):
-    """ extract locref + scmap from network """
+    """extract locref + scmap from network"""
     scmap = outputs_np[0]
     scmap = np.squeeze(scmap)
     locref = None
@@ -114,7 +114,7 @@ def multi_pose_predict(scmap, locref, stride, num_outputs):
 
 
 def getpose(image, cfg, sess, inputs, outputs, outall=False):
-    """ Extract pose """
+    """Extract pose"""
     im = np.expand_dims(image, axis=0).astype(float)
     outputs_np = sess.run(outputs, feed_dict={inputs: im})
     scmap, locref = extract_cnn_output(outputs_np, cfg)
@@ -131,7 +131,7 @@ def getpose(image, cfg, sess, inputs, outputs, outall=False):
 
 ## Functions below implement are for batch sizes > 1:
 def extract_cnn_outputmulti(outputs_np, cfg):
-    """ extract locref + scmap from network
+    """extract locref + scmap from network
     Dimensions: image batch x imagedim1 x imagedim2 x bodypart"""
     scmap = outputs_np[0]
     locref = None
@@ -163,8 +163,8 @@ def get_top_values(scmap, n_top=5):
 
 
 def getposeNP(image, cfg, sess, inputs, outputs, outall=False):
-    """ Adapted from DeeperCut, performs numpy-based faster inference on batches.
-        Introduced in https://www.biorxiv.org/content/10.1101/457242v1 """
+    """Adapted from DeeperCut, performs numpy-based faster inference on batches.
+    Introduced in https://www.biorxiv.org/content/10.1101/457242v1"""
 
     num_outputs = cfg.get("num_outputs", 1)
     outputs_np = sess.run(outputs, feed_dict={inputs: image})
diff --git a/deeplabcut/pose_estimation_tensorflow/core/predict_multianimal.py b/deeplabcut/pose_estimation_tensorflow/core/predict_multianimal.py
index 3f976c231f..b352f1ba43 100644
--- a/deeplabcut/pose_estimation_tensorflow/core/predict_multianimal.py
+++ b/deeplabcut/pose_estimation_tensorflow/core/predict_multianimal.py
@@ -18,7 +18,7 @@
 
 
 def extract_cnn_output(outputs_np, cfg):
-    """ extract locref, scmap and partaffinityfield from network """
+    """extract locref, scmap and partaffinityfield from network"""
     scmap = outputs_np[0]
     scmap = np.squeeze(scmap)
     if cfg["location_refinement"]:
@@ -39,7 +39,7 @@ def extract_cnn_output(outputs_np, cfg):
 
 
 def extract_cnn_outputmulti(outputs_np, cfg):
-    """ extract locref + scmap from network
+    """extract locref + scmap from network
     Dimensions: image batch x imagedim1 x imagedim2 x bodypart"""
     scmap = outputs_np[0]
     if cfg["location_refinement"]:
@@ -60,7 +60,13 @@ def extract_cnn_outputmulti(outputs_np, cfg):
 
 
 def compute_edge_costs(
-    pafs, peak_inds_in_batch, graph, paf_inds, n_bodyparts, n_points=10, n_decimals=3,
+    pafs,
+    peak_inds_in_batch,
+    graph,
+    paf_inds,
+    n_bodyparts,
+    n_points=10,
+    n_decimals=3,
 ):
     # Clip peak locations to PAFs dimensions
     h, w = pafs.shape[1:3]
@@ -233,7 +239,13 @@ def predict_batched_peaks_and_costs(
     )
     if peaks_gt is not None and graph:
         costs_gt = compute_edge_costs(
-            pafs, peaks_gt, graph, limbs, pose_cfg["num_joints"], n_points, n_decimals,
+            pafs,
+            peaks_gt,
+            graph,
+            limbs,
+            pose_cfg["num_joints"],
+            n_points,
+            n_decimals,
         )
         for i, costs in enumerate(costs_gt):
             preds[i]["groundtruth_costs"] = costs
@@ -268,13 +280,19 @@ def find_local_peak_indices_dilation(scmaps, radius, threshold):
     width = tf.shape(scmaps)[2]
     depth = tf.shape(scmaps)[3]
     scmaps_flat = tf.reshape(
-        tf.transpose(scmaps, [0, 3, 1, 2]), [-1, height, width, 1],
+        tf.transpose(scmaps, [0, 3, 1, 2]),
+        [-1, height, width, 1],
     )
     scmaps_dil = tf.nn.dilation2d(
-        scmaps_flat, kernel, strides=[1, 1, 1, 1], rates=[1, 1, 1, 1], padding="SAME",
+        scmaps_flat,
+        kernel,
+        strides=[1, 1, 1, 1],
+        rates=[1, 1, 1, 1],
+        padding="SAME",
     )
     scmaps_dil = tf.transpose(
-        tf.reshape(scmaps_dil, [-1, depth, height, width]), [0, 2, 3, 1],
+        tf.reshape(scmaps_dil, [-1, depth, height, width]),
+        [0, 2, 3, 1],
     )
     argmax_and_thresh_img = (scmaps > scmaps_dil) & (scmaps > threshold)
     return tf.cast(tf.where(argmax_and_thresh_img), tf.int32)
@@ -293,7 +311,10 @@ def find_local_peak_indices_skimage(scmaps, radius, threshold):
 
 
 def calc_peak_locations(
-    locrefs, peak_inds_in_batch, stride, n_decimals=3,
+    locrefs,
+    peak_inds_in_batch,
+    stride,
+    n_decimals=3,
 ):
     s, r, c, b = peak_inds_in_batch.T
     off = locrefs[s, r, c, b]
diff --git a/deeplabcut/pose_estimation_tensorflow/core/train.py b/deeplabcut/pose_estimation_tensorflow/core/train.py
index 19e8d90f34..1ad3f410da 100644
--- a/deeplabcut/pose_estimation_tensorflow/core/train.py
+++ b/deeplabcut/pose_estimation_tensorflow/core/train.py
@@ -86,7 +86,8 @@ def load_and_enqueue(sess, enqueue_op, coord, dataset, placeholders):
 def start_preloading(sess, enqueue_op, dataset, placeholders):
     coord = tf.compat.v1.train.Coordinator()
     t = threading.Thread(
-        target=load_and_enqueue, args=(sess, enqueue_op, coord, dataset, placeholders),
+        target=load_and_enqueue,
+        args=(sess, enqueue_op, coord, dataset, placeholders),
     )
     t.start()
     return coord, t
diff --git a/deeplabcut/pose_estimation_tensorflow/datasets/augmentation.py b/deeplabcut/pose_estimation_tensorflow/datasets/augmentation.py
index 7e07e271d1..5f3b9496fe 100644
--- a/deeplabcut/pose_estimation_tensorflow/datasets/augmentation.py
+++ b/deeplabcut/pose_estimation_tensorflow/datasets/augmentation.py
@@ -5,7 +5,11 @@
 
 class KeypointAwareCropToFixedSize(iaa.CropToFixedSize):
     def __init__(
-        self, width, height, max_shift=0.4, crop_sampling="hybrid",
+        self,
+        width,
+        height,
+        max_shift=0.4,
+        crop_sampling="hybrid",
     ):
         """
         Parameters
@@ -28,7 +32,9 @@ def __init__(
             or "hybrid" (alternating randomly between "uniform" and "density").
         """
         super(KeypointAwareCropToFixedSize, self).__init__(
-            width, height, name="kptscrop",
+            width,
+            height,
+            name="kptscrop",
         )
         # Clamp to 40% of crop size to ensure that at least
         # the center keypoint remains visible after the offset is applied.
diff --git a/deeplabcut/pose_estimation_tensorflow/datasets/pose_deterministic.py b/deeplabcut/pose_estimation_tensorflow/datasets/pose_deterministic.py
index 1beaa5fce3..d6175fd7e8 100644
--- a/deeplabcut/pose_estimation_tensorflow/datasets/pose_deterministic.py
+++ b/deeplabcut/pose_estimation_tensorflow/datasets/pose_deterministic.py
@@ -172,7 +172,7 @@ def make_batch(self, data_item, scale, mirror):
         im_file = data_item.im_path
         logging.debug("image %s", im_file)
         logging.debug("mirror %r", mirror)
-        image = imread(os.path.join(self.cfg["project_path"], im_file), mode="RGB")
+        image = imread(os.path.join(self.cfg["project_path"], im_file), mode="skimage")
 
         if self.has_gt:
             joints = np.copy(data_item.joints)
diff --git a/deeplabcut/pose_estimation_tensorflow/datasets/pose_imgaug.py b/deeplabcut/pose_estimation_tensorflow/datasets/pose_imgaug.py
index 5dbf1a43a7..1e61253996 100644
--- a/deeplabcut/pose_estimation_tensorflow/datasets/pose_imgaug.py
+++ b/deeplabcut/pose_estimation_tensorflow/datasets/pose_imgaug.py
@@ -300,7 +300,9 @@ def get_batch(self):
             im_file = data_item.im_path
 
             logging.debug("image %s", im_file)
-            image = imread(os.path.join(self.cfg["project_path"], im_file), mode="RGB")
+            image = imread(
+                os.path.join(self.cfg["project_path"], im_file), mode="skimage"
+            )
 
             if self.has_gt:
                 joints = np.copy(data_item.joints)
diff --git a/deeplabcut/pose_estimation_tensorflow/datasets/pose_multianimal_imgaug.py b/deeplabcut/pose_estimation_tensorflow/datasets/pose_multianimal_imgaug.py
index d6ff5f34c6..befe4f5775 100644
--- a/deeplabcut/pose_estimation_tensorflow/datasets/pose_multianimal_imgaug.py
+++ b/deeplabcut/pose_estimation_tensorflow/datasets/pose_multianimal_imgaug.py
@@ -236,7 +236,9 @@ def get_batch(self):
             im_file = data_item.im_path
 
             logging.debug("image %s", im_file)
-            image = imread(os.path.join(self.cfg["project_path"], im_file), mode="RGB")
+            image = imread(
+                os.path.join(self.cfg["project_path"], im_file), mode="skimage"
+            )
             if self.has_gt:
                 Joints = data_item.joints
                 joint_id = [
@@ -252,7 +254,12 @@ def get_batch(self):
         return batch_images, joint_ids, batch_joints, data_items
 
     def get_targetmaps_update(
-        self, joint_ids, joints, data_items, sm_size, scale,
+        self,
+        joint_ids,
+        joints,
+        data_items,
+        sm_size,
+        scale,
     ):
         part_score_targets = []
         part_score_weights = []
@@ -315,7 +322,12 @@ def calc_target_and_scoremap_sizes(self):
 
     def next_batch(self, plotting=False):
         while True:
-            (batch_images, joint_ids, batch_joints, data_items,) = self.get_batch()
+            (
+                batch_images,
+                joint_ids,
+                batch_joints,
+                data_items,
+            ) = self.get_batch()
 
             # Scale is sampled only once (per batch) to transform all of the images into same size.
             target_size, sm_size = self.calc_target_and_scoremap_sizes()
diff --git a/deeplabcut/pose_estimation_tensorflow/datasets/utils.py b/deeplabcut/pose_estimation_tensorflow/datasets/utils.py
index 1a69841782..0295eb63c2 100644
--- a/deeplabcut/pose_estimation_tensorflow/datasets/utils.py
+++ b/deeplabcut/pose_estimation_tensorflow/datasets/utils.py
@@ -31,7 +31,7 @@ def mirror_joints_map(all_joints, num_joints):
 
 
 def crop_image(joints, im, Xlabel, Ylabel, cfg):
-    """ Randomly cropping image around xlabel,ylabel taking into account size of image.
+    """Randomly cropping image around xlabel,ylabel taking into account size of image.
     Introduced in DLC 2.0 (Nature Protocols paper)"""
     widthforward = int(cfg["minsize"] + np.random.randint(cfg["rightwidth"]))
     widthback = int(cfg["minsize"] + np.random.randint(cfg["leftwidth"]))
diff --git a/deeplabcut/pose_estimation_tensorflow/export.py b/deeplabcut/pose_estimation_tensorflow/export.py
index de6f096acf..a8ff8b990e 100644
--- a/deeplabcut/pose_estimation_tensorflow/export.py
+++ b/deeplabcut/pose_estimation_tensorflow/export.py
@@ -225,7 +225,9 @@ def tf_to_pb(sess, checkpoint, output, output_dir=None):
     # create frozen graph from pbtxt file
     pb_file = os.path.normpath(output_dir + "/" + ckpt_base + ".pb")
     frozen_graph_def = tf.compat.v1.graph_util.convert_variables_to_constants(
-        sess, sess.graph_def, output,
+        sess,
+        sess.graph_def,
+        output,
     )
     with open(pb_file, "wb") as file:
         file.write(frozen_graph_def.SerializeToString())
diff --git a/deeplabcut/pose_estimation_tensorflow/lib/crossvalutils.py b/deeplabcut/pose_estimation_tensorflow/lib/crossvalutils.py
index 2cc7e1b4a6..579e885261 100644
--- a/deeplabcut/pose_estimation_tensorflow/lib/crossvalutils.py
+++ b/deeplabcut/pose_estimation_tensorflow/lib/crossvalutils.py
@@ -98,7 +98,10 @@ def _calc_separability(
 
 
 def _calc_within_between_pafs(
-    data, metadata, per_edge=True, train_set_only=True,
+    data,
+    metadata,
+    per_edge=True,
+    train_set_only=True,
 ):
     data = deepcopy(data)
     train_inds = set(metadata["data"]["trainIndices"])
@@ -319,7 +322,9 @@ def _get_n_best_paf_graphs(
         raise ValueError('`which` must be either "best" or "worst"')
 
     (within_train, _), (between_train, _) = _calc_within_between_pafs(
-        data, metadata, train_set_only=True,
+        data,
+        metadata,
+        train_set_only=True,
     )
     # Handle unlabeled bodyparts...
     existing_edges = set(k for k, v in within_train.items() if v)
@@ -397,7 +402,11 @@ def cross_validate_paf_graphs(
         cfg, params["paf_graph"]
     )
     best_graphs = _get_n_best_paf_graphs(
-        data, metadata, params["paf_graph"], ignore_inds=to_ignore, n_graphs=n_graphs,
+        data,
+        metadata,
+        params["paf_graph"],
+        ignore_inds=to_ignore,
+        n_graphs=n_graphs,
     )
     paf_scores = best_graphs[1]
     if paf_inds is None:
@@ -424,7 +433,10 @@ def cross_validate_paf_graphs(
         margin=margin,
         symmetric_kpts=symmetric_kpts,
         calibration_file=calibration_file,
-        split_inds=[metadata["data"]["trainIndices"], metadata["data"]["testIndices"],],
+        split_inds=[
+            metadata["data"]["trainIndices"],
+            metadata["data"]["testIndices"],
+        ],
     )
     # Select optimal PAF graph
     df = results[1]
diff --git a/deeplabcut/pose_estimation_tensorflow/lib/inferenceutils.py b/deeplabcut/pose_estimation_tensorflow/lib/inferenceutils.py
index 491088b381..4b09b2a9a0 100644
--- a/deeplabcut/pose_estimation_tensorflow/lib/inferenceutils.py
+++ b/deeplabcut/pose_estimation_tensorflow/lib/inferenceutils.py
@@ -841,7 +841,11 @@ def to_pickle(self, output_name):
 
 
 def calc_object_keypoint_similarity(
-    xy_pred, xy_true, sigma, margin=0, symmetric_kpts=None,
+    xy_pred,
+    xy_true,
+    sigma,
+    margin=0,
+    symmetric_kpts=None,
 ):
     visible_gt = ~np.isnan(xy_true).all(axis=1)
     if visible_gt.sum() < 2:  # At least 2 points needed to calculate scale
@@ -881,7 +885,11 @@ def calc_object_keypoint_similarity(
 
 
 def match_assemblies(
-    ass_pred, ass_true, sigma, margin=0, symmetric_kpts=None,
+    ass_pred,
+    ass_true,
+    sigma,
+    margin=0,
+    symmetric_kpts=None,
 ):
     # Only consider assemblies of at least two keypoints
     ass_pred = [a for a in ass_pred if len(a) > 1]
@@ -891,7 +899,11 @@ def match_assemblies(
     for i, a_pred in enumerate(ass_pred):
         for j, a_true in enumerate(ass_true):
             oks = calc_object_keypoint_similarity(
-                a_pred.xy, a_true.xy, sigma, margin, symmetric_kpts,
+                a_pred.xy,
+                a_true.xy,
+                sigma,
+                margin,
+                symmetric_kpts,
             )
             if ~np.isnan(oks):
                 mat[i, j] = oks
@@ -968,7 +980,11 @@ def evaluate_assembly(
         if ass_true is None:
             continue
         matched, unmatched = match_assemblies(
-            ass_pred, ass_true, oks_sigma, margin, symmetric_kpts,
+            ass_pred,
+            ass_true,
+            oks_sigma,
+            margin,
+            symmetric_kpts,
         )
         all_matched.extend(matched)
         all_unmatched.extend(unmatched)
diff --git a/deeplabcut/pose_estimation_tensorflow/lib/trackingutils.py b/deeplabcut/pose_estimation_tensorflow/lib/trackingutils.py
index 7678f1f07c..135c4a43d6 100644
--- a/deeplabcut/pose_estimation_tensorflow/lib/trackingutils.py
+++ b/deeplabcut/pose_estimation_tensorflow/lib/trackingutils.py
@@ -46,7 +46,12 @@ class BaseTracker:
     n_trackers = 0
 
     def __init__(self, dim, dim_z):
-        self.kf = kinematic_kf(dim, 1, dim_z=dim_z, order_by_dim=False,)
+        self.kf = kinematic_kf(
+            dim,
+            1,
+            dim_z=dim_z,
+            order_by_dim=False,
+        )
         self.id = self.__class__.n_trackers
         self.__class__.n_trackers += 1
         self.time_since_update = 0
diff --git a/deeplabcut/pose_estimation_tensorflow/nnets/base.py b/deeplabcut/pose_estimation_tensorflow/nnets/base.py
index cf0420f0ef..372c478e8d 100644
--- a/deeplabcut/pose_estimation_tensorflow/nnets/base.py
+++ b/deeplabcut/pose_estimation_tensorflow/nnets/base.py
@@ -76,7 +76,10 @@ def test(self, inputs):
         return self.add_inference_layers(heads)
 
     def prediction_layers(
-        self, features, scope="pose", reuse=None,
+        self,
+        features,
+        scope="pose",
+        reuse=None,
     ):
         out = {}
         n_joints = self.cfg["num_joints"]
@@ -89,26 +92,35 @@ def prediction_layers(
             )
             if self.cfg["location_refinement"]:
                 out["locref"] = prediction_layer(
-                    self.cfg, features, "locref_pred", n_joints * 2,
+                    self.cfg,
+                    features,
+                    "locref_pred",
+                    n_joints * 2,
                 )
             if (
                 self.cfg["pairwise_predict"]
                 and "multi-animal" not in self.cfg["dataset_type"]
             ):
                 out["pairwise_pred"] = prediction_layer(
-                    self.cfg, features, "pairwise_pred", n_joints * (n_joints - 1) * 2,
+                    self.cfg,
+                    features,
+                    "pairwise_pred",
+                    n_joints * (n_joints - 1) * 2,
                 )
             if (
                 self.cfg["partaffinityfield_predict"]
                 and "multi-animal" in self.cfg["dataset_type"]
             ):
                 out["pairwise_pred"] = prediction_layer(
-                    self.cfg, features, "pairwise_pred", self.cfg["num_limbs"] * 2,
+                    self.cfg,
+                    features,
+                    "pairwise_pred",
+                    self.cfg["num_limbs"] * 2,
                 )
         return out
 
     def inference(self, inputs):
-        """ Direct TF inference on GPU.
+        """Direct TF inference on GPU.
         Added with: https://arxiv.org/abs/1909.11229
         """
         heads = self.get_net(inputs)
@@ -169,7 +181,7 @@ def inference(self, inputs):
         return {"pose": pose}
 
     def add_inference_layers(self, heads):
-        """ initialized during inference """
+        """initialized during inference"""
         prob = tf.sigmoid(heads["part_pred"])
         nms_radius = int(self.cfg.get("nmsradius", 5))
 
@@ -177,16 +189,22 @@ def add_inference_layers(self, heads):
         # https://openaccess.thecvf.com/content_CVPR_2020/papers/Huang_The_Devil_Is_in_the_Details_Delving_Into_Unbiased_Data_CVPR_2020_paper.pdf
         scmaps = tf.gather(prob, tf.range(self.cfg["num_joints"]), axis=3)
         kernel = make_2d_gaussian_kernel(
-            sigma=self.cfg.get("sigma", 1), size=nms_radius * 2 + 1,
+            sigma=self.cfg.get("sigma", 1),
+            size=nms_radius * 2 + 1,
         )
         kernel = kernel[:, :, tf.newaxis, tf.newaxis]
 
         kernel_sc = tf.tile(kernel, [1, 1, tf.shape(scmaps)[3], 1])
         scmaps = tf.nn.depthwise_conv2d(
-            scmaps, kernel_sc, strides=[1, 1, 1, 1], padding="SAME",
+            scmaps,
+            kernel_sc,
+            strides=[1, 1, 1, 1],
+            padding="SAME",
         )
         peak_inds = predict_multianimal.find_local_peak_indices_maxpool_nms(
-            scmaps, nms_radius, self.cfg.get("minconfidence", 0.01),
+            scmaps,
+            nms_radius,
+            self.cfg.get("minconfidence", 0.01),
         )
         outputs = {"part_prob": prob, "peak_inds": peak_inds}
         if self.cfg["location_refinement"]:
@@ -194,7 +212,10 @@ def add_inference_layers(self, heads):
             if self.cfg.get("locref_smooth", False):
                 kernel_loc = tf.tile(kernel, [1, 1, tf.shape(locref)[3], 1])
                 locref = tf.nn.depthwise_conv2d(
-                    locref, kernel_loc, strides=[1, 1, 1, 1], padding="SAME",
+                    locref,
+                    kernel_loc,
+                    strides=[1, 1, 1, 1],
+                    padding="SAME",
                 )
             outputs["locref"] = locref
 
diff --git a/deeplabcut/pose_estimation_tensorflow/nnets/conv_blocks.py b/deeplabcut/pose_estimation_tensorflow/nnets/conv_blocks.py
index 71864b9373..4898c64cc4 100644
--- a/deeplabcut/pose_estimation_tensorflow/nnets/conv_blocks.py
+++ b/deeplabcut/pose_estimation_tensorflow/nnets/conv_blocks.py
@@ -23,19 +23,19 @@
 def _fixed_padding(inputs, kernel_size, rate=1):
     """Pads the input along the spatial dimensions independently of input size.
 
-  Pads the input such that if it was used in a convolution with 'VALID' padding,
-  the output would have the same dimensions as if the unpadded input was used
-  in a convolution with 'SAME' padding.
-
-  Args:
-    inputs: A tensor of size [batch, height_in, width_in, channels].
-    kernel_size: The kernel to be used in the conv2d or max_pool2d operation.
-    rate: An integer, rate for atrous convolution.
-
-  Returns:
-    output: A tensor of size [batch, height_out, width_out, channels] with the
-      input, either intact (if kernel_size == 1) or padded (if kernel_size > 1).
-  """
+    Pads the input such that if it was used in a convolution with 'VALID' padding,
+    the output would have the same dimensions as if the unpadded input was used
+    in a convolution with 'SAME' padding.
+
+    Args:
+      inputs: A tensor of size [batch, height_in, width_in, channels].
+      kernel_size: The kernel to be used in the conv2d or max_pool2d operation.
+      rate: An integer, rate for atrous convolution.
+
+    Returns:
+      output: A tensor of size [batch, height_out, width_out, channels] with the
+        input, either intact (if kernel_size == 1) or padded (if kernel_size > 1).
+    """
     size = kernel_size[0] + (kernel_size[0] - 1) * (rate - 1)
     kernel_size_effective = [size, size]
     pad_total = [kernel_size_effective[0] - 1, kernel_size_effective[1] - 1]
@@ -103,29 +103,29 @@ def split_separable_conv2d(
 ):
     """Separable mobilenet V1 style convolution.
 
-  Depthwise convolution, with default non-linearity,
-  followed by 1x1 depthwise convolution.  This is similar to
-  slim.separable_conv2d, but differs in tha it applies batch
-  normalization and non-linearity to depthwise. This  matches
-  the basic building of Mobilenet Paper
-  (https://arxiv.org/abs/1704.04861)
-
-  Args:
-    input_tensor: input
-    num_outputs: number of outputs
-    scope: optional name of the scope. Note if provided it will use
-    scope_depthwise for deptwhise, and scope_pointwise for pointwise.
-    normalizer_fn: which normalizer function to use for depthwise/pointwise
-    stride: stride
-    rate: output rate (also known as dilation rate)
-    endpoints: optional, if provided, will export additional tensors to it.
-    use_explicit_padding: Use 'VALID' padding for convolutions, but prepad
-      inputs so that the output dimensions are the same as if 'SAME' padding
-      were used.
-
-  Returns:
-    output tesnor
-  """
+    Depthwise convolution, with default non-linearity,
+    followed by 1x1 depthwise convolution.  This is similar to
+    slim.separable_conv2d, but differs in tha it applies batch
+    normalization and non-linearity to depthwise. This  matches
+    the basic building of Mobilenet Paper
+    (https://arxiv.org/abs/1704.04861)
+
+    Args:
+      input_tensor: input
+      num_outputs: number of outputs
+      scope: optional name of the scope. Note if provided it will use
+      scope_depthwise for deptwhise, and scope_pointwise for pointwise.
+      normalizer_fn: which normalizer function to use for depthwise/pointwise
+      stride: stride
+      rate: output rate (also known as dilation rate)
+      endpoints: optional, if provided, will export additional tensors to it.
+      use_explicit_padding: Use 'VALID' padding for convolutions, but prepad
+        inputs so that the output dimensions are the same as if 'SAME' padding
+        were used.
+
+    Returns:
+      output tesnor
+    """
 
     with _v1_compatible_scope_naming(scope) as scope:
         dw_scope = scope + "depthwise"
@@ -190,53 +190,53 @@ def expanded_conv(
 ):
     """Depthwise Convolution Block with expansion.
 
-  Builds a composite convolution that has the following structure
-  expansion (1x1) -> depthwise (kernel_size) -> projection (1x1)
-
-  Args:
-    input_tensor: input
-    num_outputs: number of outputs in the final layer.
-    expansion_size: the size of expansion, could be a constant or a callable.
-      If latter it will be provided 'num_inputs' as an input. For forward
-      compatibility it should accept arbitrary keyword arguments.
-      Default will expand the input by factor of 6.
-    stride: depthwise stride
-    rate: depthwise rate
-    kernel_size: depthwise kernel
-    residual: whether to include residual connection between input
-      and output.
-    normalizer_fn: batchnorm or otherwise
-    project_activation_fn: activation function for the project layer
-    split_projection: how many ways to split projection operator
-      (that is conv expansion->bottleneck)
-    split_expansion: how many ways to split expansion op
-      (that is conv bottleneck->expansion) ops will keep depth divisible
-      by this value.
-    split_divisible_by: make sure every split group is divisible by this number.
-    expansion_transform: Optional function that takes expansion
-      as a single input and returns output.
-    depthwise_location: where to put depthwise covnvolutions supported
-      values None, 'input', 'output', 'expansion'
-    depthwise_channel_multiplier: depthwise channel multiplier:
-    each input will replicated (with different filters)
-    that many times. So if input had c channels,
-    output will have c x depthwise_channel_multpilier.
-    endpoints: An optional dictionary into which intermediate endpoints are
-      placed. The keys "expansion_output", "depthwise_output",
-      "projection_output" and "expansion_transform" are always populated, even
-      if the corresponding functions are not invoked.
-    use_explicit_padding: Use 'VALID' padding for convolutions, but prepad
-      inputs so that the output dimensions are the same as if 'SAME' padding
-      were used.
-    padding: Padding type to use if `use_explicit_padding` is not set.
-    scope: optional scope.
-
-  Returns:
-    Tensor of depth num_outputs
-
-  Raises:
-    TypeError: on inval
-  """
+    Builds a composite convolution that has the following structure
+    expansion (1x1) -> depthwise (kernel_size) -> projection (1x1)
+
+    Args:
+      input_tensor: input
+      num_outputs: number of outputs in the final layer.
+      expansion_size: the size of expansion, could be a constant or a callable.
+        If latter it will be provided 'num_inputs' as an input. For forward
+        compatibility it should accept arbitrary keyword arguments.
+        Default will expand the input by factor of 6.
+      stride: depthwise stride
+      rate: depthwise rate
+      kernel_size: depthwise kernel
+      residual: whether to include residual connection between input
+        and output.
+      normalizer_fn: batchnorm or otherwise
+      project_activation_fn: activation function for the project layer
+      split_projection: how many ways to split projection operator
+        (that is conv expansion->bottleneck)
+      split_expansion: how many ways to split expansion op
+        (that is conv bottleneck->expansion) ops will keep depth divisible
+        by this value.
+      split_divisible_by: make sure every split group is divisible by this number.
+      expansion_transform: Optional function that takes expansion
+        as a single input and returns output.
+      depthwise_location: where to put depthwise covnvolutions supported
+        values None, 'input', 'output', 'expansion'
+      depthwise_channel_multiplier: depthwise channel multiplier:
+      each input will replicated (with different filters)
+      that many times. So if input had c channels,
+      output will have c x depthwise_channel_multpilier.
+      endpoints: An optional dictionary into which intermediate endpoints are
+        placed. The keys "expansion_output", "depthwise_output",
+        "projection_output" and "expansion_transform" are always populated, even
+        if the corresponding functions are not invoked.
+      use_explicit_padding: Use 'VALID' padding for convolutions, but prepad
+        inputs so that the output dimensions are the same as if 'SAME' padding
+        were used.
+      padding: Padding type to use if `use_explicit_padding` is not set.
+      scope: optional scope.
+
+    Returns:
+      Tensor of depth num_outputs
+
+    Raises:
+      TypeError: on inval
+    """
     with tf.compat.v1.variable_scope(
         scope, default_name="expanded_conv"
     ) as s, tf.compat.v1.name_scope(s.original_name_scope):
@@ -339,20 +339,20 @@ def expanded_conv(
 def split_conv(input_tensor, num_outputs, num_ways, scope, divisible_by=8, **kwargs):
     """Creates a split convolution.
 
-  Split convolution splits the input and output into
-  'num_blocks' blocks of approximately the same size each,
-  and only connects $i$-th input to $i$ output.
-
-  Args:
-    input_tensor: input tensor
-    num_outputs: number of output filters
-    num_ways: num blocks to split by.
-    scope: scope for all the operators.
-    divisible_by: make sure that every part is divisiable by this.
-    **kwargs: will be passed directly into conv2d operator
-  Returns:
-    tensor
-  """
+    Split convolution splits the input and output into
+    'num_blocks' blocks of approximately the same size each,
+    and only connects $i$-th input to $i$ output.
+
+    Args:
+      input_tensor: input tensor
+      num_outputs: number of output filters
+      num_ways: num blocks to split by.
+      scope: scope for all the operators.
+      divisible_by: make sure that every part is divisiable by this.
+      **kwargs: will be passed directly into conv2d operator
+    Returns:
+      tensor
+    """
     b = input_tensor.get_shape().as_list()[3]
 
     if num_ways == 1 or min(b // num_ways, num_outputs // num_ways) < divisible_by:
diff --git a/deeplabcut/pose_estimation_tensorflow/nnets/layers.py b/deeplabcut/pose_estimation_tensorflow/nnets/layers.py
index feff5544ac..ee7f3056df 100644
--- a/deeplabcut/pose_estimation_tensorflow/nnets/layers.py
+++ b/deeplabcut/pose_estimation_tensorflow/nnets/layers.py
@@ -56,5 +56,10 @@ def prediction_layer_stage(cfg, input, name, num_outputs):
         weights_regularizer=slim.l2_regularizer(cfg["weight_decay"]),
     ):
         with tf.compat.v1.variable_scope(name):
-            pred = slim.conv2d(input, num_outputs, kernel_size=[3, 3], stride=1,)
+            pred = slim.conv2d(
+                input,
+                num_outputs,
+                kernel_size=[3, 3],
+                stride=1,
+            )
             return pred
diff --git a/deeplabcut/pose_estimation_tensorflow/nnets/mobilenet.py b/deeplabcut/pose_estimation_tensorflow/nnets/mobilenet.py
index 87952bd8de..1c18b3686d 100644
--- a/deeplabcut/pose_estimation_tensorflow/nnets/mobilenet.py
+++ b/deeplabcut/pose_estimation_tensorflow/nnets/mobilenet.py
@@ -69,9 +69,17 @@ def extract_features(self, inputs):
         return net, end_points
 
     def prediction_layers(
-        self, features, end_points, scope="pose", reuse=None,
+        self,
+        features,
+        end_points,
+        scope="pose",
+        reuse=None,
     ):
-        out = super(PoseMobileNet, self).prediction_layers(features, scope, reuse,)
+        out = super(PoseMobileNet, self).prediction_layers(
+            features,
+            scope,
+            reuse,
+        )
         with tf.compat.v1.variable_scope(scope, reuse=reuse):
             if self.cfg["intermediate_supervision"]:
                 out["part_pred_interm"] = prediction_layer(
diff --git a/deeplabcut/pose_estimation_tensorflow/nnets/multi.py b/deeplabcut/pose_estimation_tensorflow/nnets/multi.py
index 092c7d265d..81358ea942 100644
--- a/deeplabcut/pose_estimation_tensorflow/nnets/multi.py
+++ b/deeplabcut/pose_estimation_tensorflow/nnets/multi.py
@@ -97,7 +97,10 @@ def prediction_layer(cfg, input, name, num_outputs):
     ):
         with tf.compat.v1.variable_scope(name):
             pred = slim.conv2d_transpose(
-                input, num_outputs, kernel_size=[3, 3], stride=2,
+                input,
+                num_outputs,
+                kernel_size=[3, 3],
+                stride=2,
             )
             return pred
 
@@ -141,7 +144,12 @@ def extract_features(self, inputs):
         return net, end_points
 
     def prediction_layers(
-        self, features, end_points, input_shape, scope="pose", reuse=None,
+        self,
+        features,
+        end_points,
+        input_shape,
+        scope="pose",
+        reuse=None,
     ):
         net_type = self.cfg["net_type"]
         if self.cfg["multi_stage"]:  # MuNet! (multi_stage decoder + multi_fusion)
@@ -399,7 +407,11 @@ def prediction_layers(
                         scope="block4",
                     )
             net = tf.concat([bank_3, upsampled_features], 3)
-            out = super(PoseMultiNet, self).prediction_layers(net, scope, reuse,)
+            out = super(PoseMultiNet, self).prediction_layers(
+                net,
+                scope,
+                reuse,
+            )
             with tf.compat.v1.variable_scope(scope, reuse=reuse):
                 if (
                     self.cfg["intermediate_supervision"]
diff --git a/deeplabcut/pose_estimation_tensorflow/nnets/resnet.py b/deeplabcut/pose_estimation_tensorflow/nnets/resnet.py
index 90fb716bc2..fcced2dc59 100644
--- a/deeplabcut/pose_estimation_tensorflow/nnets/resnet.py
+++ b/deeplabcut/pose_estimation_tensorflow/nnets/resnet.py
@@ -37,14 +37,25 @@ def extract_features(self, inputs):
         im_centered = self.center_inputs(inputs)
         with slim.arg_scope(resnet_v1.resnet_arg_scope()):
             net, end_points = net_fun(
-                im_centered, global_pool=False, output_stride=16, is_training=False,
+                im_centered,
+                global_pool=False,
+                output_stride=16,
+                is_training=False,
             )
         return net, end_points
 
     def prediction_layers(
-        self, features, end_points, scope="pose", reuse=None,
+        self,
+        features,
+        end_points,
+        scope="pose",
+        reuse=None,
     ):
-        out = super(PoseResnet, self).prediction_layers(features, scope, reuse,)
+        out = super(PoseResnet, self).prediction_layers(
+            features,
+            scope,
+            reuse,
+        )
         with tf.compat.v1.variable_scope(scope, reuse=reuse):
             if self.cfg["intermediate_supervision"]:
                 layer_name = "resnet_v1_{}/block{}/unit_{}/bottleneck_v1"
diff --git a/deeplabcut/pose_estimation_tensorflow/predict_multianimal.py b/deeplabcut/pose_estimation_tensorflow/predict_multianimal.py
index b75bf816a0..eaa66e3047 100644
--- a/deeplabcut/pose_estimation_tensorflow/predict_multianimal.py
+++ b/deeplabcut/pose_estimation_tensorflow/predict_multianimal.py
@@ -36,7 +36,7 @@ def AnalyzeMultiAnimalVideo(
     robust_nframes=False,
     use_shelve=False,
 ):
-    """ Helper function for analyzing a video with multiple individuals """
+    """Helper function for analyzing a video with multiple individuals"""
 
     print("Starting to analyze % ", video)
     vname = Path(video).stem
@@ -99,7 +99,14 @@ def AnalyzeMultiAnimalVideo(
             )
         else:
             PredicteData, nframes = GetPoseandCostsS(
-                cfg, dlc_cfg, sess, inputs, outputs, vid, nframes, shelf_path,
+                cfg,
+                dlc_cfg,
+                sess,
+                inputs,
+                outputs,
+                vid,
+                nframes,
+                shelf_path,
             )
 
         stop = time.time()
@@ -138,9 +145,17 @@ def AnalyzeMultiAnimalVideo(
 
 
 def GetPoseandCostsF(
-    cfg, dlc_cfg, sess, inputs, outputs, cap, nframes, batchsize, shelf_path,
+    cfg,
+    dlc_cfg,
+    sess,
+    inputs,
+    outputs,
+    cap,
+    nframes,
+    batchsize,
+    shelf_path,
 ):
-    """ Batchwise prediction of pose """
+    """Batchwise prediction of pose"""
     strwidth = int(np.ceil(np.log10(nframes)))  # width for strings
     batch_ind = 0  # keeps track of which image within a batch should be written to
     batch_num = 0  # keeps track of which batch you are at
@@ -156,7 +171,10 @@ def GetPoseandCostsF(
     inds = []
 
     if shelf_path:
-        db = shelve.open(shelf_path, protocol=pickle.DEFAULT_PROTOCOL,)
+        db = shelve.open(
+            shelf_path,
+            protocol=pickle.DEFAULT_PROTOCOL,
+        )
     else:
         db = dict()
     db["metadata"] = {
@@ -187,7 +205,11 @@ def GetPoseandCostsF(
             inds.append(counter)
             if batch_ind == batchsize - 1:
                 D = predict.predict_batched_peaks_and_costs(
-                    dlc_cfg, frames, sess, inputs, outputs,
+                    dlc_cfg,
+                    frames,
+                    sess,
+                    inputs,
+                    outputs,
                 )
                 for ind, data in zip(inds, D):
                     db["frame" + str(ind).zfill(strwidth)] = data
@@ -200,7 +222,11 @@ def GetPoseandCostsF(
         elif counter >= nframes:
             if batch_ind > 0:
                 D = predict.predict_batched_peaks_and_costs(
-                    dlc_cfg, frames, sess, inputs, outputs,
+                    dlc_cfg,
+                    frames,
+                    sess,
+                    inputs,
+                    outputs,
                 )
                 for ind, data in zip(inds, D):
                     db["frame" + str(ind).zfill(strwidth)] = data
@@ -219,13 +245,16 @@ def GetPoseandCostsF(
 
 
 def GetPoseandCostsS(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes, shelf_path):
-    """ Non batch wise pose estimation for video cap."""
+    """Non batch wise pose estimation for video cap."""
     strwidth = int(np.ceil(np.log10(nframes)))  # width for strings
     if cfg["cropping"]:
         cap.set_bbox(cfg["x1"], cfg["x2"], cfg["y1"], cfg["y2"])
 
     if shelf_path:
-        db = shelve.open(shelf_path, protocol=pickle.DEFAULT_PROTOCOL,)
+        db = shelve.open(
+            shelf_path,
+            protocol=pickle.DEFAULT_PROTOCOL,
+        )
     else:
         db = dict()
     db["metadata"] = {
@@ -255,7 +284,11 @@ def GetPoseandCostsS(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes, shelf_pa
             if frame.shape[-1] == 4:
                 frame = rgba2rgb(frame)
             dets = predict.predict_batched_peaks_and_costs(
-                dlc_cfg, np.expand_dims(frame, axis=0), sess, inputs, outputs,
+                dlc_cfg,
+                np.expand_dims(frame, axis=0),
+                sess,
+                inputs,
+                outputs,
             )
             db[key] = dets[0]
             del dets
diff --git a/deeplabcut/pose_estimation_tensorflow/predict_videos.py b/deeplabcut/pose_estimation_tensorflow/predict_videos.py
index 3e828ca33e..bc21581c03 100644
--- a/deeplabcut/pose_estimation_tensorflow/predict_videos.py
+++ b/deeplabcut/pose_estimation_tensorflow/predict_videos.py
@@ -57,11 +57,11 @@ def analyze_videos(
     modelprefix="",
     robust_nframes=False,
     allow_growth=False,
+    use_shelve=False,
     auto_track=True,
     n_tracks=None,
     calibrate=False,
     identity_only=False,
-    use_shelve=False,
 ):
     """
     Makes prediction based on a trained network. The index of the trained network is specified by parameters in the config file (in particular the variable 'snapshotindex')
@@ -126,14 +126,21 @@ def analyze_videos(
     allow_growth: bool, default false.
         For some smaller GPUs the memory issues happen. If true, the memory allocator does not pre-allocate the entire specified
         GPU memory region, instead starting small and growing as needed. See issue: https://forum.image.sc/t/how-to-stop-running-out-of-vram/30551/2
+    
+    use_shelve: bool, optional (default=False)
+        By default, data are dumped in a pickle file at the end of the video analysis.
+        Otherwise, data are written to disk on the fly using a "shelf"; i.e., a pickle-based,
+        persistent, database-like object by default, resulting in constant memory footprint.
 
+    The following parameters are only relevant for multi-animal projects:
+    
     auto_track: bool, optional (default=True)
         By default, tracking and stitching are automatically performed, producing the final h5 data file.
-        This is equivalent to the behavior of single-animal projects.
+        This is equivalent to the behavior for single-animal projects.
 
         If False, one must run `convert_detections2tracklets` and `stitch_tracklets` afterwards, in order to obtain the h5 file.
 
-        This function has 3 related sub-calls:
+    This function has 3 related sub-calls:
 
     identity_only: bool, optional (default=False)
         If True and animal identity was learned by the model,
@@ -150,12 +157,6 @@ def analyze_videos(
         passed if the number of animals in the video is different from
         the number of animals the model was trained on.
 
-
-    use_shelve: bool, optional (default=False)
-        By default, data are dumped in a pickle file at the end of the video analysis.
-        Otherwise, data are written to disk on the fly using a "shelf"; i.e., a pickle-based,
-        persistent, database-like object by default, resulting in constant memory footprint.
-
     Examples
     --------
 
@@ -427,7 +428,7 @@ def checkcropping(cfg, cap):
 
 
 def GetPoseF(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes, batchsize):
-    """ Batchwise prediction of pose """
+    """Batchwise prediction of pose"""
     PredictedData = np.zeros(
         (nframes, dlc_cfg["num_outputs"] * 3 * len(dlc_cfg["all_joints_names"]))
     )
@@ -479,7 +480,7 @@ def GetPoseF(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes, batchsize):
 
 
 def GetPoseS(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes):
-    """ Non batch wise pose estimation for video cap."""
+    """Non batch wise pose estimation for video cap."""
     if cfg["cropping"]:
         ny, nx = checkcropping(cfg, cap)
 
@@ -517,7 +518,7 @@ def GetPoseS(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes):
 
 
 def GetPoseS_GTF(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes):
-    """ Non batch wise pose estimation for video cap."""
+    """Non batch wise pose estimation for video cap."""
     if cfg["cropping"]:
         ny, nx = checkcropping(cfg, cap)
 
@@ -562,7 +563,7 @@ def GetPoseS_GTF(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes):
 
 
 def GetPoseF_GTF(cfg, dlc_cfg, sess, inputs, outputs, cap, nframes, batchsize):
-    """ Batchwise prediction of pose """
+    """Batchwise prediction of pose"""
     PredictedData = np.zeros((nframes, 3 * len(dlc_cfg["all_joints_names"])))
     batch_ind = 0  # keeps track of which image within a batch should be written to
     batch_num = 0  # keeps track of which batch you are at
@@ -633,7 +634,7 @@ def getboundingbox(x, y, nx, ny, margin):
 def GetPoseDynamic(
     cfg, dlc_cfg, sess, inputs, outputs, cap, nframes, detectiontreshold, margin
 ):
-    """ Non batch wise pose estimation for video cap by dynamically cropping around previously detected parts."""
+    """Non batch wise pose estimation for video cap by dynamically cropping around previously detected parts."""
     if cfg["cropping"]:
         ny, nx = checkcropping(cfg, cap)
     else:
@@ -717,7 +718,7 @@ def AnalyzeVideo(
     TFGPUinference=True,
     dynamic=(False, 0.5, 10),
 ):
-    """ Helper function for analyzing a video. """
+    """Helper function for analyzing a video."""
     print("Starting to analyze % ", video)
 
     if destfolder is None:
@@ -845,17 +846,13 @@ def AnalyzeVideo(
 
 
 def GetPosesofFrames(
-    cfg, dlc_cfg, sess, inputs, outputs, directory, framelist, nframes, batchsize, rgb
+    cfg, dlc_cfg, sess, inputs, outputs, directory, framelist, nframes, batchsize
 ):
-    """ Batchwise prediction of pose for frame list in directory"""
-    # from skimage.io import imread
+    """Batchwise prediction of pose for frame list in directory"""
     from deeplabcut.utils.auxfun_videos import imread
 
     print("Starting to extract posture")
-    if rgb:
-        im = imread(os.path.join(directory, framelist[0]), mode="RGB")
-    else:
-        im = imread(os.path.join(directory, framelist[0]))
+    im = imread(os.path.join(directory, framelist[0]), mode="skimage")
 
     ny, nx, nc = np.shape(im)
     print(
@@ -897,11 +894,7 @@ def GetPosesofFrames(
 
     if batchsize == 1:
         for counter, framename in enumerate(framelist):
-            # frame=imread(os.path.join(directory,framename),mode='RGB')
-            if rgb:
-                im = imread(os.path.join(directory, framename), mode="RGB")
-            else:
-                im = imread(os.path.join(directory, framename))
+            im = imread(os.path.join(directory, framename), mode="skimage")
 
             if counter % step == 0:
                 pbar.update(step)
@@ -920,10 +913,7 @@ def GetPosesofFrames(
             (batchsize, ny, nx, 3), dtype="ubyte"
         )  # this keeps all the frames of a batch
         for counter, framename in enumerate(framelist):
-            if rgb:
-                im = imread(os.path.join(directory, framename), mode="RGB")
-            else:
-                im = imread(os.path.join(directory, framename))
+            im = imread(os.path.join(directory, framename), mode="skimage")
 
             if counter % step == 0:
                 pbar.update(step)
@@ -967,7 +957,6 @@ def analyze_time_lapse_frames(
     trainingsetindex=0,
     gputouse=None,
     save_as_csv=False,
-    rgb=True,
     modelprefix="",
 ):
     """
@@ -1003,9 +992,6 @@ def analyze_time_lapse_frames(
     save_as_csv: bool, optional
         Saves the predictions in a .csv file. The default is ``False``; if provided it must be either ``True`` or ``False``
 
-    rbg: bool, optional.
-        Whether to load image as rgb; Note e.g. some tiffs do not alow that option in imread, then just set this to false.
-
     Examples
     --------
     If you want to analyze all frames in /analysis/project/timelapseexperiment1
@@ -1142,7 +1128,6 @@ def analyze_time_lapse_frames(
                     framelist,
                     nframes,
                     dlc_cfg["batch_size"],
-                    rgb,
                 )
                 stop = time.time()
 
@@ -1189,7 +1174,13 @@ def analyze_time_lapse_frames(
 
 
 def _convert_detections_to_tracklets(
-    cfg, inference_cfg, data, metadata, output_path, greedy=False, calibrate=False,
+    cfg,
+    inference_cfg,
+    data,
+    metadata,
+    output_path,
+    greedy=False,
+    calibrate=False,
 ):
     track_method = cfg.get("default_track_method", "ellipse")
     if track_method not in ("box", "skeleton", "ellipse"):
diff --git a/deeplabcut/pose_estimation_tensorflow/training.py b/deeplabcut/pose_estimation_tensorflow/training.py
index 4b0906e395..a135f03092 100644
--- a/deeplabcut/pose_estimation_tensorflow/training.py
+++ b/deeplabcut/pose_estimation_tensorflow/training.py
@@ -13,7 +13,7 @@
 
 
 def return_train_network_path(config, shuffle=1, trainingsetindex=0, modelprefix=""):
-    """ Returns the training and test pose config file names as well as the folder where the snapshot is
+    """Returns the training and test pose config file names as well as the folder where the snapshot is
     Parameters
     ----------
     config : string
diff --git a/deeplabcut/pose_estimation_tensorflow/visualizemaps.py b/deeplabcut/pose_estimation_tensorflow/visualizemaps.py
index be5f0aa3f0..b51e7c69e7 100644
--- a/deeplabcut/pose_estimation_tensorflow/visualizemaps.py
+++ b/deeplabcut/pose_estimation_tensorflow/visualizemaps.py
@@ -226,14 +226,14 @@ def extract_maps(
             DATA = {}
             for imageindex, imagename in tqdm(Indices):
                 image = imread(
-                    os.path.join(cfg["project_path"], *imagename), mode="RGB"
+                    os.path.join(cfg["project_path"], *imagename), mode="skimage"
                 )
-                if image.shape[-1] == 4:
-                    image = rgba2rgb(image)
+
                 if scale != 1:
                     image = imresize(image, scale)
 
                 image_batch = data_to_input(image)
+
                 # Compute prediction with the CNN
                 outputs_np = sess.run(outputs, feed_dict={inputs: image_batch})
 
diff --git a/deeplabcut/refine_training_dataset/outlier_frames.py b/deeplabcut/refine_training_dataset/outlier_frames.py
index 9079753b9f..3ccc1f508c 100644
--- a/deeplabcut/refine_training_dataset/outlier_frames.py
+++ b/deeplabcut/refine_training_dataset/outlier_frames.py
@@ -428,8 +428,8 @@ def extract_outlier_frames(
 
 
 def convertparms2start(pn):
-    """ Creating a start value for sarimax in case of an value error
-    See: https://groups.google.com/forum/#!topic/pystatsmodels/S_Fo53F25Rk """
+    """Creating a start value for sarimax in case of an value error
+    See: https://groups.google.com/forum/#!topic/pystatsmodels/S_Fo53F25Rk"""
     if "ar." in pn:
         return 0
     elif "ma." in pn:
@@ -484,8 +484,8 @@ def FitSARIMAXModel(x, p, pcutoff, alpha, ARdegree, MAdegree, nforecast=0, disp=
 def compute_deviations(
     Dataframe, dataname, p_bound, alpha, ARdegree, MAdegree, storeoutput=None
 ):
-    """ Fits Seasonal AutoRegressive Integrated Moving Average with eXogenous regressors model to data and computes confidence interval
-    as well as mean fit. """
+    """Fits Seasonal AutoRegressive Integrated Moving Average with eXogenous regressors model to data and computes confidence interval
+    as well as mean fit."""
 
     print("Fitting state-space models with parameters:", ARdegree, MAdegree)
     df_x, df_y, df_likelihood = Dataframe.values.reshape((Dataframe.shape[0], -1, 3)).T
@@ -790,7 +790,7 @@ def PlottingSingleFrame(
     strwidth=4,
     savelabeled=True,
 ):
-    """ Label frame and save under imagename / this is already cropped (for clip) """
+    """Label frame and save under imagename / this is already cropped (for clip)"""
     from skimage import io
 
     imagename1 = os.path.join(tmpfolder, "img" + str(index).zfill(strwidth) + ".png")
@@ -859,7 +859,7 @@ def PlottingSingleFramecv2(
     strwidth=4,
     savelabeled=True,
 ):
-    """ Label frame and save under imagename / cap is not already cropped. """
+    """Label frame and save under imagename / cap is not already cropped."""
     from skimage import io
 
     imagename1 = os.path.join(tmpfolder, "img" + str(index).zfill(strwidth) + ".png")
diff --git a/deeplabcut/refine_training_dataset/stitch.py b/deeplabcut/refine_training_dataset/stitch.py
index b2820d1425..1d503fd7b0 100644
--- a/deeplabcut/refine_training_dataset/stitch.py
+++ b/deeplabcut/refine_training_dataset/stitch.py
@@ -573,7 +573,10 @@ def compute_max_gap(tracklets):
         return max_gap
 
     def build_graph(
-        self, nodes=None, max_gap=None, weight_func=None,
+        self,
+        nodes=None,
+        max_gap=None,
+        weight_func=None,
     ):
         if nodes is None:
             nodes = self.tracklets
diff --git a/deeplabcut/utils/auxfun_videos.py b/deeplabcut/utils/auxfun_videos.py
index 400a13c006..98e49131b3 100644
--- a/deeplabcut/utils/auxfun_videos.py
+++ b/deeplabcut/utils/auxfun_videos.py
@@ -8,13 +8,20 @@
 https://github.com/AlexEMG/DeepLabCut/blob/master/AUTHORS
 Licensed under GNU Lesser General Public License v3.0
 """
+
+#from deeplabcut.utils.auxfun_videos import imread
+#auxfun_videos.imread(image_path, mode="skimage")
+
+import skimage.color
+from skimage import io
+from skimage.util import img_as_ubyte
 import cv2
 import datetime
 import numpy as np
 import os
 import subprocess
-import warnings
-
+import warnings                        
+                            
 
 class VideoReader:
     def __init__(self, video_path):
@@ -345,10 +352,20 @@ def check_video_integrity(video_path):
     vid.check_integrity()
     vid.check_integrity_robust()
 
-
-# Historically DLC used: from scipy.misc import imread, imresize >> deprecated functions
-def imread(path, mode=None):
-    return cv2.imread(path, cv2.IMREAD_UNCHANGED)[..., ::-1]  # ~10% faster than using cv2.cvtColor
+def imread(image_path, mode="skimage"):
+    ''' Read image either with skimage or cv2. 
+    Returns frame in uint with 3 color channels. '''
+    if mode == "skimage":
+        image = io.imread(image_path)
+        if image.ndim == 2 or image.shape[-1] == 1:
+            image = skimage.color.gray2rgb(image)
+        elif image.shape[-1] == 4:
+            image = skimage.color.rgba2rgb(image)
+
+        return img_as_ubyte(image)
+    
+    elif mode=="cv2":
+        return cv2.imread(image_path, cv2.IMREAD_UNCHANGED)[..., ::-1]  # ~10% faster than using cv2.cvtColor
 
 
 # https://docs.opencv.org/3.4.0/da/d54/group__imgproc__transform.html#ga5bb5a1fea74ea38e1a5445ca803ff121
diff --git a/deeplabcut/utils/skeleton.py b/deeplabcut/utils/skeleton.py
index 678baea2a3..e01d92fcfb 100644
--- a/deeplabcut/utils/skeleton.py
+++ b/deeplabcut/utils/skeleton.py
@@ -75,6 +75,7 @@ def __init__(self, config_path):
         sep = "/" if "/" in row else "\\"
         if sep != os.path.sep:
             row = row.replace(sep, os.path.sep)
+        
         self.image = io.imread(os.path.join(self.cfg["project_path"], row))
         self.inds = set()
         self.segs = set()
diff --git a/deeplabcut/version.py b/deeplabcut/version.py
index 5e6da8d866..61a3c219be 100644
--- a/deeplabcut/version.py
+++ b/deeplabcut/version.py
@@ -8,5 +8,5 @@
 Licensed under GNU Lesser General Public License v3.0
 """
 
-__version__ = "2.2.0.3"
+__version__ = "2.2.0.4"
 VERSION = __version__
diff --git a/examples/test.sh b/examples/test.sh
index 29bb214261..8b54f90f63 100755
--- a/examples/test.sh
+++ b/examples/test.sh
@@ -6,7 +6,7 @@ rm -r OUT
 cd ..
 pip uninstall deeplabcut
 python3 setup.py sdist bdist_wheel
-pip install dist/deeplabcut-2.2.0.3-py3-none-any.whl
+pip install dist/deeplabcut-2.2.0.4-py3-none-any.whl
 
 cd examples
 
diff --git a/examples/testM1.sh b/examples/testM1.sh
index c3df2174ae..9c2efca044 100755
--- a/examples/testM1.sh
+++ b/examples/testM1.sh
@@ -6,7 +6,7 @@ rm -r OUT
 cd ..
 pip uninstall deeplabcut
 pythonw setup.py sdist bdist_wheel
-pip install dist/deeplabcut-2.2.0.3-py3-none-any.whl
+pip install dist/deeplabcut-2.2.0.4-py3-none-any.whl
 
 # download: https://drive.google.com/file/d/17pSwfoNuyf3YR8vCaVggHeI-pMQ3xL7l/view?usp=sharing
 # assuming it's in Downloads...
diff --git a/examples/testscript.py b/examples/testscript.py
index adfc39934d..34e4aa5524 100644
--- a/examples/testscript.py
+++ b/examples/testscript.py
@@ -129,7 +129,10 @@
     # Check the training image paths are correctly stored as arrays of strings
     trainingsetfolder = auxiliaryfunctions.GetTrainingSetFolder(cfg)
     datafile, _ = auxiliaryfunctions.GetDataandMetaDataFilenames(
-        trainingsetfolder, 0.8, 1, cfg,
+        trainingsetfolder,
+        0.8,
+        1,
+        cfg,
     )
     mlab = sio.loadmat(os.path.join(cfg["project_path"], datafile))["dataset"]
     num_images = mlab.shape[1]
@@ -366,7 +369,10 @@ def make_frame(t):
     print("will be used for 3D testscript...")
     # TENSORPACK could fail in WINDOWS...
     deeplabcut.create_training_dataset(
-        path_config_file, Shuffles=[2], net_type=NET, augmenter_type=augmenter_type3,
+        path_config_file,
+        Shuffles=[2],
+        net_type=NET,
+        augmenter_type=augmenter_type3,
     )
 
     posefile = os.path.join(
diff --git a/examples/testscript_multianimal.py b/examples/testscript_multianimal.py
index fc96bd5b04..5df60776b6 100644
--- a/examples/testscript_multianimal.py
+++ b/examples/testscript_multianimal.py
@@ -111,7 +111,10 @@
     # Check the training image paths are correctly stored as arrays of strings
     trainingsetfolder = auxiliaryfunctions.GetTrainingSetFolder(cfg)
     datafile, _ = auxiliaryfunctions.GetDataandMetaDataFilenames(
-        trainingsetfolder, 0.8, 1, cfg,
+        trainingsetfolder,
+        0.8,
+        1,
+        cfg,
     )
     datafile = datafile.split(".mat")[0] + ".pickle"
     with open(os.path.join(cfg["project_path"], datafile), "rb") as f:
@@ -178,7 +181,9 @@
 
     print("Convert detections to tracklets...")
     deeplabcut.convert_detections2tracklets(
-        config_path, [new_video_path], "mp4",
+        config_path,
+        [new_video_path],
+        "mp4",
     )
     print("Tracklets created...")
 
@@ -191,25 +196,36 @@
 
     print("Plotting trajectories...")
     deeplabcut.plot_trajectories(
-        config_path, [new_video_path], "mp4",
+        config_path,
+        [new_video_path],
+        "mp4",
     )
     print("Trajectory plotted.")
 
     print("Creating labeled video...")
     deeplabcut.create_labeled_video(
-        config_path, [new_video_path], "mp4", save_frames=False, color_by="individual",
+        config_path,
+        [new_video_path],
+        "mp4",
+        save_frames=False,
+        color_by="individual",
     )
     print("Labeled video created.")
 
     print("Filtering predictions...")
     deeplabcut.filterpredictions(
-        config_path, [new_video_path], "mp4",
+        config_path,
+        [new_video_path],
+        "mp4",
     )
     print("Predictions filtered.")
 
     print("Extracting outlier frames...")
     deeplabcut.extract_outlier_frames(
-        config_path, [new_video_path], "mp4", automatic=True,
+        config_path,
+        [new_video_path],
+        "mp4",
+        automatic=True,
     )
     print("Outlier frames extracted.")
 
diff --git a/examples/testscript_openfielddata_augmentationcomparison.py b/examples/testscript_openfielddata_augmentationcomparison.py
index ebc5161667..88c33834d3 100644
--- a/examples/testscript_openfielddata_augmentationcomparison.py
+++ b/examples/testscript_openfielddata_augmentationcomparison.py
@@ -2,151 +2,10 @@
 # -*- coding: utf-8 -*-
 """
 
-This is a test script to compare the loaders. tensorpack allows much more choices for augmentation. The parameters
-can be set in pose_dataset_tensorpack.py and of course specifically in each pose_config.yaml file before training. In fact,
-pose_dataset_tensorpack.py will fall back to default parameters if they are not defined in pose_config.yaml and one is
-using dataset_type:'tensorpack'
+This is a test script to compare the loaders and models. 
 
-This script creates one identical split for the openfield test dataset and trains it with the
-standard loader and the tensorpack loader for k iterations in DLC 2.0 docker with TF 1.8 on a NVIDIA GTX 1080Ti.
-
-My results were (Run with DLC 2.0.9 in Sept 2019)
-
-**With standard loader:**
-
-Training iterations:	%Training dataset	Shuffle number	 Train error(px)	 Test error(px)	p-cutoff used	Train error with p-cutoff	Test error with p-cutoff
-10000	80	2	2.64	3.11	0.4	2.64	3.11
-20000	80	2	2.26	2.72	0.4	2.26	2.72
-30000	80	2	1.71	2.28	0.4	1.71	2.28
-40000	80	2	1.88	2.61	0.4	1.88	2.61
-50000	80	2	1.86	2.32	0.4	1.86	2.32
-60000	80	2	1.92	2.42	0.4	1.92	2.42
-70000	80	2	2.38	3.04	0.4	2.38	3.04
-80000	80	2	1.55	2.34	0.4	1.55	2.34
-90000	80	2	1.5	2.27	0.4	1.5	2.27
-100000	80	2	1.52	2.34	0.4	1.52	2.34
-
-
-**With tensorpack loader:**
-
-Training iterations:	%Training dataset	Shuffle number	 Train error(px)	 Test error(px)	p-cutoff used	Train error with p-cutoff	Test error with p-cutoff
-10000	80	3	2.35	2.91	0.4	2.35	2.91
-20000	80	3	3.28	3.51	0.4	3.28	3.51
-30000	80	3	1.57	2.24	0.4	1.57	2.24
-40000	80	3	3.54	4.17	0.4	3.54	4.17
-50000	80	3	1.76	2.74	0.4	1.76	2.74
-60000	80	3	2.85	3.39	0.4	2.85	3.39
-70000	80	3	3.88	4.71	0.4	3.88	4.71
-80000	80	3	1.2	2.06	0.4	1.2	2.06
-90000	80	3	2.2	3.07	0.4	2.2	3.07
-100000	80	3	1.06	1.96	0.4	1.06	1.96
-
-
-For details on TensorPack check out:
-
-A Neural Net Training Interface on TensorFlow, with focus on speed + flexibility
-https://github.com/tensorpack/tensorpack
-
-My results were (Run with DLC 2.2b5 in May 2020) for 20k iterations
-
-Imagaug augmentation:
-
-Results for 20000  training iterations: 95 1 train error: 3.25 pixels. Test error: 4.98  pixels.
-With pcutoff of 0.4  train error: 3.25 pixels. Test error: 4.98 pixels
-
-Default augmentation:
-
-Results for 20000  training iterations: 95 2 train error: 2.5 pixels. Test error: 4.08  pixels.
-With pcutoff of 0.4  train error: 2.5 pixels. Test error: 4.08 pixels
-
-Tensorpack augmentation:
-
-Results for 20000  training iterations: 95 3 train error: 3.06 pixels. Test error: 4.78  pixels.
-With pcutoff of 0.4  train error: 3.06 pixels. Test error: 4.78 pixels
-
-My results were (Run with DLC *2.2b7* in July 2020) for 20k iterations
-
-Attention: default changed!
-
-***Default = Imagaug**** augmentation:
-
-Done and results stored for snapshot:  snapshot-20000
-Results for 20000  training iterations: 95 1 train error: 2.93 pixels. Test error: 3.09  pixels.
-With pcutoff of 0.4  train error: 2.93 pixels. Test error: 3.09 pixels
-
-Scalecrop (was = default) augmentation:
-
-Done and results stored for snapshot:  snapshot-20000
-Results for 20000  training iterations: 95 2 train error: 2.5 pixels. Test error: 2.57  pixels.
-With pcutoff of 0.4  train error: 2.5 pixels. Test error: 2.57 pixels
-
-Tensorpack augmentation:
-
-Done and results stored for snapshot:  snapshot-20000
-Results for 20000  training iterations: 95 3 train error: 3.1 pixels. Test error: 3.29  pixels.
-With pcutoff of 0.4  train error: 3.1 pixels. Test error: 3.29 pixels
-
-My results were (Run with DLC *2.2b7* on August 1st 2020) for 10k iterations
-
-Imgaug:
-Results for 10000  training iterations: 95 1 train error: 3.78 pixels. Test error: 3.89  pixels.
-With pcutoff of 0.4  train error: 3.78 pixels. Test error: 3.89 pixels
-
-Scalecrop:
-Done and results stored for snapshot:  snapshot-10000
-Results for 10000  training iterations: 95 2 train error: 2.81 pixels. Test error: 2.46  pixels.
-With pcutoff of 0.4  train error: 2.81 pixels. Test error: 2.46 pixels
-
-Tensorpack:
-Done and results stored for snapshot:  snapshot-10000
-Results for 10000  training iterations: 95 3 train error: 3.76 pixels. Test error: 3.98  pixels.
-With pcutoff of 0.4  train error: 3.76 pixels. Test error: 3.98 pixels
-
-
-My results were (Run with DLC *2.2b8* on Sept 7 2020) for 10k iterations
-
-Imgaug:
-Results for 10000  training iterations: 95 1 train error: 2.63 pixels. Test error: 3.88  pixels.
-With pcutoff of 0.4  train error: 2.63 pixels. Test error: 3.88 pixels
-
-Scalecrop:
-Results for 10000  training iterations: 95 2 train error: 3.08 pixels. Test error: 4.02  pixels.
-With pcutoff of 0.4  train error: 3.08 pixels. Test error: 4.02 pixels
-
-Tensorpack:
-Results for 10000  training iterations: 95 3 train error: 2.9 pixels. Test error: 3.31  pixels.
-With pcutoff of 0.4  train error: 2.9 pixels. Test error: 3.31 pixels
-
-My results were (Run with DLC *2.1.9* in Jan 2021) for 10 k iterations
-
-**ResNet50
-Imgaug:
-Results for 100000  training iterations: 95 1 train error: 2.13 pixels. Test error: 2.22  pixels.
-With pcutoff of 0.4  train error: 2.13 pixels. Test error: 2.22 pixels
-
-Scalecrop:
-Results for 100000  training iterations: 95 2 train error: 1.47 pixels. Test error: 1.77  pixels.
-With pcutoff of 0.4  train error: 1.47 pixels. Test error: 1.77 pixels
-
-Tensorpack:
-Results for 100000  training iterations: 95 3 train error: 2.09 pixels. Test error: 2.36  pixels.
-With pcutoff of 0.4  train error: 2.09 pixels. Test error: 2.36 pixels
-
-**EffNet-b3
-Imgaug:
-Results for 100000  training iterations: 95 4 train error: 2.39 pixels. Test error: 2.57  pixels.
-With pcutoff of 0.4  train error: 2.39 pixels. Test error: 2.57 pixels
-
-Scalecrop:
-Results for 100000  training iterations: 95 5 train error: 2.26 pixels. Test error: 2.24  pixels.
-With pcutoff of 0.4  train error: 2.26 pixels. Test error: 2.24 pixels
-
-Tensorpack:
-Results for 100000  training iterations: 95 6 train error: 1.65 pixels. Test error: 2.24  pixels.
-With pcutoff of 0.4  train error: 1.65 pixels. Test error: 2.24 pixels
-
-Notice: despite the higher RMSE for imgaug due to the augmentation,
-the network performs much better on the testvideo (see Neuron Primer: https://www.cell.com/neuron/pdf/S0896-6273(20)30717-0.pdf)
+This script creates one identical splits for the openfield test dataset and trains it with imgaug (default), scalecrop
+and the tensorpack loader. We also compare 3 backbones (mobilenet, resnet, efficientnet)
 
 My results were (Run with DLC *2.10.4* in Apr 2021) for 100 k iterations
 
@@ -165,6 +24,25 @@
 Results for 100000  training iterations: 95 3 train error: 1.35 pixels. Test error: 2.3  pixels.
 With pcutoff of 0.4  train error: 1.35 pixels. Test error: 2.3 pixels
 
+Results Jan 2020:
+
+MobileNetV2 0.35
+Results for 10000  training iterations: 95 1 train error: 5.79 pixels. Test error: 5.63  pixels.
+With pcutoff of 0.4  train error: 5.79 pixels. Test error: 5.63 pixels
+
+ResNet 50
+Results for 10000  training iterations: 95 2 train error: 3.61 pixels. Test error: 3.7  pixels.
+With pcutoff of 0.4  train error: 3.61 pixels. Test error: 3.7 pixels
+
+EffNet-b3
+Results for 10000  training iterations: 95 3 train error: 6.86 pixels. Test error: 6.63  pixels.
+With pcutoff of 0.4  train error: 6.86 pixels. Test error: 6.63 pixels
+
+TODO: Note we should still optimize the MobNet & EffNet learning rates for this dataset (also training is pretty short!)
+
+Notice: despite the higher RMSE for imgaug due to the augmentation,
+the network performs much better on the testvideo (see Neuron Primer: https://www.cell.com/neuron/pdf/S0896-6273(20)30717-0.pdf)
+
 """
 
 
@@ -183,34 +61,39 @@
 displayiters = 500
 Shuffles = 1 + np.arange(6)
 
-deeplabcut.load_demo_data(path_config_file)
-## Create one split and make Shuffle 2 and 3 have the same split.
+deeplabcut.load_demo_data(path_config_file, createtrainingset=False)
+## Create one identical splits for 3 networks and 3 augmentations
+
 ###Note that the new function in DLC 2.1 simplifies network/augmentation comparisons greatly:
-deeplabcut.create_training_model_comparison(
+Shuffles = deeplabcut.create_training_model_comparison(
     path_config_file,
     num_shuffles=1,
-    net_types=["resnet_50", "efficientnet-b3"],
+    net_types=["mobilenet_v2_0.35", "resnet_50", "efficientnet-b3"],
     augmenter_types=["imgaug", "scalecrop", "tensorpack"],
 )
 
-
-for shuffle in Shuffles:
-
+for idx, shuffle in enumerate(Shuffles):
     posefile, _, _ = deeplabcut.return_train_network_path(
         path_config_file, shuffle=shuffle
     )
 
-    edits = {"decay_steps": maxiters, "lr_init": 0.0005}  # * 8}  # for EfficientNet
-    DLC_config = deeplabcut.auxiliaryfunctions.edit_config(posefile, edits)
-
-    if shuffle % 3 == 1:  # imgaug
+    # Setting specific parameters for training
+    if idx % 3 == 0:  # imgaug
         edits = {"rotation": 180, "motion_blur": True}
         DLC_config = deeplabcut.auxiliaryfunctions.edit_config(posefile, edits)
-
-    elif shuffle % 3 == 0:  # Tensorpack:
+    elif idx % 3 == 2:  # Tensorpack
         edits = {"rotation": 180, "noise_sigma": 0.01}
         DLC_config = deeplabcut.auxiliaryfunctions.edit_config(posefile, edits)
 
+    if idx > 5:  # EfficientNet
+        print(posefile, "changing now!!")
+        edits = {
+            "decay_steps": maxiters,
+            "lr_init": 0.0005,
+        }
+        DLC_config = deeplabcut.auxiliaryfunctions.edit_config(posefile, edits)
+
+for shuffle in Shuffles:
     print("TRAIN NETWORK", shuffle)
     deeplabcut.train_network(
         path_config_file,
@@ -221,10 +104,6 @@
         max_snapshots_to_keep=11,
     )
 
-print("EVALUATE")
-deeplabcut.evaluate_network(path_config_file, Shuffles=Shuffles, plotting=True)
-
-for shuffles in Shuffle:
     print("Analyze Video")
 
     videofile_path = os.path.join(
@@ -236,3 +115,6 @@
     print("Create Labeled Video and plot")
     deeplabcut.create_labeled_video(path_config_file, [videofile_path], shuffle=shuffle)
     deeplabcut.plot_trajectories(path_config_file, [videofile_path], shuffle=shuffle)
+
+print("EVALUATE")
+deeplabcut.evaluate_network(path_config_file, Shuffles=Shuffles, plotting=False)
diff --git a/examples/testscript_openfielddata_netcomparison.py b/examples/testscript_openfielddata_netcomparison.py
deleted file mode 100644
index e09353f249..0000000000
--- a/examples/testscript_openfielddata_netcomparison.py
+++ /dev/null
@@ -1,89 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-This is a test script to compare the networks. On Jan 3rd 2020:
-
-Jan 2020:
-MobileNetV2 0.35
-Results for 10000  training iterations: 95 1 train error: 5.79 pixels. Test error: 5.63  pixels.
-With pcutoff of 0.4  train error: 5.79 pixels. Test error: 5.63 pixels
-
-ResNet 50
-Results for 10000  training iterations: 95 2 train error: 3.61 pixels. Test error: 3.7  pixels.
-With pcutoff of 0.4  train error: 3.61 pixels. Test error: 3.7 pixels
-
-EffNet-b3
-Results for 10000  training iterations: 95 3 train error: 6.86 pixels. Test error: 6.63  pixels.
-With pcutoff of 0.4  train error: 6.86 pixels. Test error: 6.63 pixels
-
-Note: Not too good on video either!
-
-TODO: Note we should still optimize the MobNet & EffNet learning rates for this dataset (also training is pretty short!)
-TODO: change to frozen backbone!
-"""
-
-
-import os
-
-os.environ["CUDA_VISIBLE_DEVICES"] = str(0)
-import deeplabcut
-import numpy as np
-
-# Loading example data set
-path_config_file = os.path.join(os.getcwd(), "openfield-Pranav-2018-10-30/config.yaml")
-cfg = deeplabcut.auxiliaryfunctions.read_config(path_config_file)
-maxiters = 10000
-
-deeplabcut.load_demo_data(path_config_file)
-
-## Create one split and make Shuffle 2 and 3 have the same split.
-###Note that the new function in DLC 2.1 simplifies network/augmentation comparisons greatly:
-deeplabcut.create_training_model_comparison(
-    path_config_file,
-    num_shuffles=1,
-    net_types=["mobilenet_v2_0.35", "resnet_50", "efficientnet-b3"],
-    augmenter_types=["imgaug"],
-)
-
-freezeencoder = False  # True
-for shuffle in 1 + np.arange(3):
-
-    posefile, _, _ = deeplabcut.return_train_network_path(
-        path_config_file, shuffle=shuffle
-    )
-
-    # for EfficientNet
-    edits = {
-        "decay_steps": maxiters,
-        "lr_init": 0.0005 * 12,
-        "freezeencoder": freezeencoder,
-    }
-    DLC_config = deeplabcut.auxiliaryfunctions.edit_config(posefile, edits)
-    # imgaug
-    edits = {"rotation": 180, "motion_blur": True, "freezeencoder": freezeencoder}
-    DLC_config = deeplabcut.auxiliaryfunctions.edit_config(posefile, edits)
-
-    print("TRAIN NETWORK", shuffle)
-    deeplabcut.train_network(
-        path_config_file,
-        shuffle=shuffle,
-        saveiters=10000,
-        displayiters=200,
-        maxiters=maxiters,
-        max_snapshots_to_keep=11,
-    )
-
-    print("EVALUATE")
-    deeplabcut.evaluate_network(path_config_file, Shuffles=[shuffle], plotting=True)
-
-    print("Analyze Video")
-
-    videofile_path = os.path.join(
-        os.getcwd(), "openfield-Pranav-2018-10-30", "videos", "m3v1mp4.mp4"
-    )
-
-    deeplabcut.analyze_videos(path_config_file, [videofile_path], shuffle=shuffle)
-
-    print("Create Labeled Video and plot")
-    deeplabcut.create_labeled_video(path_config_file, [videofile_path], shuffle=shuffle)
-    deeplabcut.plot_trajectories(path_config_file, [videofile_path], shuffle=shuffle)
diff --git a/reinstall.sh b/reinstall.sh
index 24a46d5a2e..c01280aef8 100755
--- a/reinstall.sh
+++ b/reinstall.sh
@@ -1,3 +1,3 @@
 pip uninstall deeplabcut
 python3 setup.py sdist bdist_wheel
-pip install dist/deeplabcut-2.2.0.3-py3-none-any.whl
+pip install dist/deeplabcut-2.2.0.4-py3-none-any.whl
diff --git a/setup.py b/setup.py
index eb01d2f134..a4f219b30f 100644
--- a/setup.py
+++ b/setup.py
@@ -16,7 +16,7 @@
 
 setuptools.setup(
     name="deeplabcut",
-    version="2.2.0.3",
+    version="2.2.0.4",
     author="A. & M. Mathis Labs",
     author_email="alexander@deeplabcut.org",
     description="Markerless pose-estimation of user-defined features with deep learning",
diff --git a/tests/test_trainingsetmanipulation.py b/tests/test_trainingsetmanipulation.py
index 2e9578ff65..99b0de0bf3 100644
--- a/tests/test_trainingsetmanipulation.py
+++ b/tests/test_trainingsetmanipulation.py
@@ -10,13 +10,15 @@
     trainingsetmanipulation,
     multiple_individuals_trainingsetmanipulation,
 )
+
+from deeplabcut.utils.auxfun_videos import imread
 from deeplabcut.utils.conversioncode import guarantee_multiindex_rows
-from skimage import io, color
+from skimage import color, io
 
 
 def test_read_image_shape_fast(tmp_path):
     path_rgb_image = os.path.join(TEST_DATA_DIR, "image.png")
-    img = io.imread(path_rgb_image)
+    img = imread(path_rgb_image, mode="skimage")
     shape = img.shape
     assert read_image_shape_fast(path_rgb_image) == (shape[2], shape[0], shape[1])
     path_gray_image = str(tmp_path / "gray.png")