docs: def set_data(self, mjd, mag, err, whitten=False, **kwarg):
          """
          Saves the light curve data, where 
          mjd: modified julian date (time instants)
          mag: stellar magnitude 
          err: photometric error
          
          If whitten is True the data is normalized to have zero
          mean and unit standard deviation. Note that robust 
          estimators of these quantities are used
  
          TODO: verify that arrays are non empty, non constant, etc
          """
          
          self.mjd = mjd
          self.N = len(mjd)
          weights = np.power(err, -2.0)
          self.weights = weights/np.sum(weights)
          self.scale = robust_scale(mag, self.weights)
          if whitten:
              self.mag = (mag - robust_center(mag, self.weights))/self.scale
              self.err = err/self.scale
          else:
              self.mag = mag
              self.err = err
          self.mjd = self.mjd.astype('float32')
          self.mag = self.mag.astype('float32')
          self.err = self.err.astype('float32')
          if self.method == 'QMICS' or self.method == 'QMIEU' or self.method == 'QME':
              if whitten:
                  hm = 0.9*self.N**(-0.2) # Silverman's rule, assuming data is whittened
              else:
                  hm = 0.9*self.scale*self.N**(-0.2)
              if 'h_KDE_M' in kwarg:
                  hm = hm*kwarg['h_KDE_M']
              hp = 1.0 # How to choose this more appropietly?
              if 'h_KDE_P' in kwarg:
                  hp = hp*kwarg['h_KDE_P']
              kernel = 0  # Select the kernel for the magnitudes, 0 is safe
              if 'kernel' in kwarg:
                  kernel = kwarg['kernel']
              if self.debug:
                  print("Kernel bandwidths: %f , %f" %(hm, hp))
              self.my_QMI = QMI(self.mjd, self.mag, self.err, hm, hp, kernel)
          elif self.method == 'LKSL':  # Lafler-Kinman Minimum String Length
              self.my_SL = LKSL(self.mjd, self.mag, self.err)
          elif self.method == 'PDM1':  # Phase Dispersion Minimization
              Nbins = int(self.N/3)
              if 'Nbins' in kwarg:
                  Nbins = kwarg['Nbins']
              self.my_PDM = PDM(self.mjd, self.mag, self.err, Nbins)
          elif self.method == 'MHAOV':  # Orthogonal Multiharmonics AoV periodogram
              Nharmonics = 1
              if 'Nharmonics' in kwarg:
                  Nharmonics = kwarg["Nharmonics"]
              self.my_AOV = AOV(self.mjd, self.mag, self.err, Nharmonics)