(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.
update internal variables for sampling the distribution with the
        current covariance matrix C. This method is O(N^3), if C is not diagonal.
update internal variables for sampling the distribution with the current covariance matrix C . This method is O ( N^3 ) if C is not diagonal .
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.sp.neg.cmuexp:  # cave:
            self.update_exponential(self.Zneg, -self.sp.neg.cmuexp)
            # self.C += self.Zpos  # pos update after Zneg would be the correct update, overall:
            # self.C = self.Zpos + Cs * Mh.expms(-self.sp.neg.cmuexp*Csi*self.Zneg*Csi) * Cs
            self.Zneg = np.zeros((self.N, self.N))

        if self.sigma_vec is not 1 and not np.all(self.sigma_vec == 1):
            self.C = dot(dot(np.diag(self.sigma_vec), self.C), np.diag(self.sigma_vec))
            self.sigma_vec[:] = 1

        if self.opts['CMA_const_trace'] in (True, 1, 2):  # normalize trace of C
            if self.opts['CMA_const_trace'] == 2:
                s = np.exp(np.mean(np.log(self.dC)))
            else:
                s = np.mean(self.dC)
            self.C /= s
            self.dC /= s
        self.C = (self.C + self.C.T) / 2
        # 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

        if type(self.opts['CMA_eigenmethod']) == type(1):
            print('WARNING: option CMA_eigenmethod should be a function, not an integer')
            if self.opts['CMA_eigenmethod'] == -1:
                # pygsl
                # easy to install (well, in Windows install gsl binaries first,
                # set system path to respective libgsl-0.dll (or cp the dll to
                # python\DLLS ?), in unzipped pygsl edit
                # gsl_dist/gsl_site_example.py into gsl_dist/gsl_site.py
                # and run "python setup.py build" and "python setup.py install"
                # in MINGW32)
                if 1 < 3:  # import pygsl on the fly
                    try:
                        import pygsl.eigen.eigenvectors  # TODO efficient enough?
                    except ImportError:
                        print('WARNING: could not find pygsl.eigen module, either install pygsl \n' +
                              '  or set option CMA_eigenmethod=1 (is much slower), option set to 1')
                        self.opts['CMA_eigenmethod'] = 0  # use 0 if 1 is too slow

                    self.D, self.B = pygsl.eigen.eigenvectors(self.C)

            elif self.opts['CMA_eigenmethod'] == 0:
                # TODO: thoroughly test np.linalg.eigh
                #       numpy.linalg.eig crashes in 200-D
                #       and EVecs with same EVals are not orthogonal
                self.D, self.B = np.linalg.eigh(self.C)  # self.B[i] is a row and not an eigenvector
            else:  # is overall two;ten times slower in 10;20-D
                self.D, self.B = Misc.eig(self.C)  # def eig, see below
        else:
            self.D, self.B = self.opts['CMA_eigenmethod'](self.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)
        idx = np.argsort(self.D)
        self.D = self.D[idx]
        self.B = self.B[:,idx]  # self.B[i] is a row, columns self.B[:,i] are eigenvectors
        # assert(all(self.B[self.countiter % self.N] == self.B[self.countiter % self.N,:]))

        # qqqqqqqqqq
        if 11 < 3:  # limit condition number to 1e13
            climit = 1e13  # cave: conditioncov termination is 1e14
            if self.D[-1] / self.D[0] > climit:
                self.D += self.D[-1] / climit
            for i in xrange(self.N):
                self.C[i][i] += self.D[-1] / climit

        if 11 < 3 and any(abs(sum(self.B[:,0:self.N-1] * self.B[:,1:], 0)) > 1e-6):
            print('B is not orthogonal')
            print(self.D)
            print(sum(self.B[:,0:self.N-1] * self.B[:,1:], 0))
        else:
            # is O(N^3)
            # assert(sum(abs(self.C - np.dot(self.D * self.B,  self.B.T))) < N**2*1e-11)
            pass
        self.D **= 0.5
        self.itereigenupdated = self.countiter
def updateBD ( self ) : # itereigenupdated is always up-to-date in the diagonal case # just double check here if self . itereigenupdated == self . countiter : return if self . sp . neg . cmuexp : # cave: self . update_exponential ( self . Zneg , - self . sp . neg . cmuexp ) # self.C += self.Zpos  # pos update after Zneg would be the correct update, overall: # self.C = self.Zpos + Cs * Mh.expms(-self.sp.neg.cmuexp*Csi*self.Zneg*Csi) * Cs self . Zneg = np . zeros ( ( self . N , self . N ) ) if self . sigma_vec is not 1 and not np . all ( self . sigma_vec == 1 ) : self . C = dot ( dot ( np . diag ( self . sigma_vec ) , self . C ) , np . diag ( self . sigma_vec ) ) self . sigma_vec [ : ] = 1 if self . opts [ 'CMA_const_trace' ] in ( True , 1 , 2 ) : # normalize trace of C if self . opts [ 'CMA_const_trace' ] == 2 : s = np . exp ( np . mean ( np . log ( self . dC ) ) ) else : s = np . mean ( self . dC ) self . C /= s self . dC /= s self . C = ( self . C + self . C . T ) / 2 # 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 if type ( self . opts [ 'CMA_eigenmethod' ] ) == type ( 1 ) : print ( 'WARNING: option CMA_eigenmethod should be a function, not an integer' ) if self . opts [ 'CMA_eigenmethod' ] == - 1 : # pygsl # easy to install (well, in Windows install gsl binaries first, # set system path to respective libgsl-0.dll (or cp the dll to # python\DLLS ?), in unzipped pygsl edit # gsl_dist/gsl_site_example.py into gsl_dist/gsl_site.py # and run "python setup.py build" and "python setup.py install" # in MINGW32) if 1 < 3 : # import pygsl on the fly try : import pygsl . eigen . eigenvectors # TODO efficient enough? except ImportError : print ( 'WARNING: could not find pygsl.eigen module, either install pygsl \n' + '  or set option CMA_eigenmethod=1 (is much slower), option set to 1' ) self . opts [ 'CMA_eigenmethod' ] = 0 # use 0 if 1 is too slow self . D , self . B = pygsl . eigen . eigenvectors ( self . C ) elif self . opts [ 'CMA_eigenmethod' ] == 0 : # TODO: thoroughly test np.linalg.eigh #       numpy.linalg.eig crashes in 200-D #       and EVecs with same EVals are not orthogonal self . D , self . B = np . linalg . eigh ( self . C ) # self.B[i] is a row and not an eigenvector else : # is overall two;ten times slower in 10;20-D self . D , self . B = Misc . eig ( self . C ) # def eig, see below else : self . D , self . B = self . opts [ 'CMA_eigenmethod' ] ( self . 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) idx = np . argsort ( self . D ) self . D = self . D [ idx ] self . B = self . B [ : , idx ] # self.B[i] is a row, columns self.B[:,i] are eigenvectors # assert(all(self.B[self.countiter % self.N] == self.B[self.countiter % self.N,:])) # qqqqqqqqqq if 11 < 3 : # limit condition number to 1e13 climit = 1e13 # cave: conditioncov termination is 1e14 if self . D [ - 1 ] / self . D [ 0 ] > climit : self . D += self . D [ - 1 ] / climit for i in xrange ( self . N ) : self . C [ i ] [ i ] += self . D [ - 1 ] / climit if 11 < 3 and any ( abs ( sum ( self . B [ : , 0 : self . N - 1 ] * self . B [ : , 1 : ] , 0 ) ) > 1e-6 ) : print ( 'B is not orthogonal' ) print ( self . D ) print ( sum ( self . B [ : , 0 : self . N - 1 ] * self . B [ : , 1 : ] , 0 ) ) else : # is O(N^3) # assert(sum(abs(self.C - np.dot(self.D * self.B,  self.B.T))) < N**2*1e-11) pass self . D **= 0.5 self . itereigenupdated = self . countiter