docs: def updateBD(self):
          """update internal variables for sampling the distribution with the
          current covariance matrix C. This method is O(N^3), if C is not diagonal.
  
          """
          # itereigenupdated is always up-to-date in the diagonal case
          # just double check here
          if self.itereigenupdated == self.countiter:
              return
          if self.opts['CMA_diagonal'] >= self.countiter:
              _print_warning("updateBD called in CMA_diagonal mode, " +
                             "this should be considered a bug", "updateBD",
                             iteration=self.countiter)
          # C has already positive updates, here come the additional negative updates
          if self.sp.neg.cmuexp:
              C_shrunken = (1 - self.sp.cmu - self.sp.c1)**(self.countiter - self.itereigenupdated)
              clip_fac = 0.60  # 0.9 is sufficient to prevent degeneration in small dimension
              if hasattr(self.opts['vv'], '__getitem__') and self.opts['vv'][0] == 'sweep_ccov_neg':
                  clip_fac = 0.98
              if (self.countiter - self.itereigenupdated) * self.sp.neg.cmuexp * self.N \
                  < clip_fac * C_shrunken:
                  # pos.def. guarantied, because vectors are normalized
                  self.C -= self.sp.neg.cmuexp * self._Yneg
              else:
                  max_warns = 1
                  try:
                      self._updateBD_warnings += 1
                  except AttributeError:
                      self._updateBD_warnings = 1
                  if self.opts['verbose'] > 1 and \
                                  self._updateBD_warnings <= max_warns:
                      _print_warning('doing two additional eigen' +
                                     'decompositions to guarantee pos.def.',
                                     'updateBD', 'CMAEvolutionStrategy')
                      if self._updateBD_warnings == max_warns:
                          _print_warning('further warnings are surpressed',
                                         'updateBD')
                  self.decompose_C()
                  _tmp_inverse_root_C = dot(self.B / self.D, self.B.T)
                  _tmp_inverse_root_C = (_tmp_inverse_root_C + _tmp_inverse_root_C.T) / 2
                  Zneg = dot(dot(_tmp_inverse_root_C, self._Yneg), _tmp_inverse_root_C)
                  eigvals, eigvecs = self.opts['CMA_eigenmethod'](Zneg)
                  self.count_eigen += 1
                  if max(eigvals) * self.sp.neg.cmuexp <= clip_fac:
                     self.C -= self.sp.neg.cmuexp * self._Yneg
                  elif 1 < 3:
                      self.C -= (clip_fac / max(eigvals)) * self._Yneg
                      _print_warning(
                          'clipped learning rate for negative weights, ' +
                          'maximal eigenvalue = %f, maxeig * ccov = %f > %f'
                           % (max(eigvals), max(eigvals) * self.sp.neg.cmuexp, clip_fac),
                          iteration=self.countiter)
                      if 1 < 3: # let's check
                          eigvals, eigvecs = self.opts['CMA_eigenmethod'](self.C)
                          self.count_eigen += 1
                          print('new min eigenval = %e, old = %e'
                                % (min(eigvals), min(self.D)**2))
                          if min(eigvals) > 0:
                              print('new cond = %e, old = %e'
                                    % (max(eigvals) / min(eigvals),
                                       (max(self.D) / min(self.D))**2))
                  else: # guaranties pos.def. unconditionally
                      _print_warning('exponential update for negative weights (internally more expensive)',
                                     iteration=self.countiter)
                      self.update_exponential(self._Yneg, -self.sp.neg.cmuexp)
                      # self.C = self.Ypos + Cs * Mh.expms(-self.sp.neg.cmuexp*Csi*self.Yneg*Csi) * Cs
              # Yneg = self.Yneg  # for temporary debugging, can be removed
              self._Yneg = np.zeros((self.N, self.N))
  
          if hasattr(self.opts['vv'], '__getitem__') and self.opts['vv'][0].startswith('sweep_ccov'):
              self.opts['CMA_const_trace'] = True
          if self.opts['CMA_const_trace'] in (True, 1, 2):  # normalize trace of C
              if self.opts['CMA_const_trace'] == 2:
                  s = np.exp(2 * np.mean(np.log(self.D)))  # or geom average of dC?
              else:
                  s = np.mean(np.diag(self.C))
              self.C /= s
  
          dC = np.diag(self.C)
          if max(dC) / min(dC) > 1e8:
              # allows for much larger condition numbers, if axis-parallel
              self.sigma_vec *= np.diag(self.C)**0.5
              self.C = self.correlation_matrix()
              _print_warning('condition in coordinate system exceeded 1e8' +
                  ', rescaled to 1')
  
          # self.C = np.triu(self.C) + np.triu(self.C,1).T  # should work as well
          # self.D, self.B = eigh(self.C) # hermitian, ie symmetric C is assumed
  
          self.decompose_C()
          # assert(sum(self.D-DD) < 1e-6)
          # assert(sum(sum(np.dot(BB, BB.T)-np.eye(self.N))) < 1e-6)
          # assert(sum(sum(np.dot(BB * DD, BB.T) - self.C)) < 1e-6)
          # assert(all(self.B[self.countiter % self.N] == self.B[self.countiter % self.N,:]))
  
          # qqqqqqqqqq
          # is O(N^3)
          # assert(sum(abs(self.C - np.dot(self.D * self.B,  self.B.T))) < N**2*1e-11)
  
          if 1 < 3 and max(self.D) / min(self.D) > 1e6 and self.gp.isidentity:
              # TODO: allow to do this again
              # dmean_prev = dot(self.B, (1. / self.D) * dot(self.B.T, (self.mean - 0*self.mean_old) / self.sigma_vec))
              self.gp._tf_matrix = (self.sigma_vec * dot(self.B * self.D, self.B.T).T).T
              self.gp._tf_matrix_inv = (dot(self.B / self.D, self.B.T).T / self.sigma_vec).T
              self.gp.tf_pheno = lambda x: dot(self.gp._tf_matrix, x)
              self.gp.tf_geno = lambda x: dot(self.gp._tf_matrix_inv, x)  # not really necessary
              self.gp.isidentity = False
              assert self.mean is not self.mean_old
              self.mean = self.gp.geno(self.mean)  # same as tf_geno
              self.mean_old = self.gp.geno(self.mean_old)  # not needed?
              self.pc = self.gp.geno(self.pc)
              self.D[:] = 1.0
              self.B = np.eye(self.N)
              self.C = np.eye(self.N)
              self.dC[:] = 1.0
              self.sigma_vec = 1
              # dmean_now = dot(self.B, (1. / self.D) * dot(self.B.T, (self.mean - 0*self.mean_old) / self.sigma_vec))
              # assert Mh.vequals_approximately(dmean_now, dmean_prev)
              _print_warning('\n geno-pheno transformation introduced based on current C,\n injected solutions become "invalid" in this iteration',
                             'updateBD', 'CMAEvolutionStrategy', self.countiter)
  
          self.itereigenupdated = self.countiter