maf.py 文件源码

def processing(self):
        '''
        Perform the maf transformation
        '''
        # covariance of original bands
        sigma = numpy.ma.cov(self.variates_stack.T, allow_masked=True)  # TODO: Jenkins can't run this lines with allow_masked =True, whyyyy??
        # covariance for horizontal and vertical shifts
        sigmadh = numpy.ma.cov(self.H.T, allow_masked=True)
        sigmadv = numpy.ma.cov(self.V.T, allow_masked=True)

        sigma = numpy.cov(self.variates_stack.T)
        # covariance for horizontal and vertical shifts
        sigmadh = numpy.cov(self.H.T)
        sigmadv = numpy.cov(self.V.T)      

        # simple pooling of shifts
        sigmad = 0.5 * (numpy.array(sigmadh) + numpy.array(sigmadv))
        # evalues, vec1 = scipy.linalg.eig(sigmad, sigma)
        evalues, vec1 = linalg.eig(sigmad, sigma)

        # Sort eigen values from smallest to largest and apply this order to
        # eigen vector matrix
        sort_index = evalues.argsort()  
        evalues = evalues[sort_index]
        vec1 = vec1[:, sort_index]
        # autocorrelation
        # ac= 1-0.5*vec1
        HH = 1 / numpy.std(self.variates_stack, 0, ddof=1)
        diagvariates = numpy.diag(HH)
        invstderrmaf = numpy.diag((1 / numpy.sqrt(numpy.diag(vec1.T * sigma * vec1))))
        HHH = numpy.zeros((self.number_of_maf_variates), float)
        for b in range(self.number_of_maf_variates):  
            # logger.info("Calculating component %d of MAF transformation" % b)  
            HHH[b] = cmp(numpy.sum((diagvariates * sigma * vec1 * invstderrmaf)[b]), 0)
        sgn = numpy.diag(HHH)  # assure positiviy
        v = numpy.dot(vec1, sgn)
        N = numpy.shape(self.variates_stack)[0]
        X = self.variates_stack - numpy.tile(numpy.mean(self.variates_stack, 0), (N, 1))
        # scale v to give MAFs with unit variance
        aux1 = numpy.dot(numpy.dot(v.T, sigma), v)  # dispersion of MAFs
        aux2 = 1 / numpy.sqrt(numpy.diag(aux1))
        aux3 = numpy.tile(aux2.T, (self.number_of_maf_variates, 1))
        v = v * aux3  # now dispersion is unit matrix
        self.mafs = numpy.dot(X, v)