code: def layer_covariance(layer1, layer2=None):
      
      layer2 = layer2 or layer1
      act1, act2 = layer1.activations, layer2.activations
      num_datapoints = act1.shape[0]  # cast to avoid numpy type promotion during division
      return np.matmul(act1.T, act2) / float(num_datapoints)
docs: Computes the covariance matrix between the neurons of two layers. If only one
      layer is passed, computes the symmetric covariance matrix of that layer.