code: def mad(a, axis=None, c=1.4826, return_med=False):
      
      a = checkma(a)
      #return np.ma.median(np.fabs(a - np.ma.median(a))) / c
      if a.count() > 0:
          if axis is None:
              med = fast_median(a)
              out = fast_median(np.fabs(a - med)) * c
          else:
              med = np.ma.median(a, axis=axis)
              #This is necessary for broadcasting
              med = np.expand_dims(med, axis=axis)
              out = np.ma.median(np.ma.fabs(a - med), axis=axis) * c
      else:
          out = np.ma.masked
      if return_med:
          out = (out, med)
      return out
docs: Compute normalized median absolute difference
     
      Can also return median array, as this can be expensive, and often we want both med and nmad
  
      Note: 1.4826 = 1/0.6745