Python scipy.signal.correlate2d() Examples

def _erase_motif_occurrences(self, seqs_onehot, ppm, ppm_bg, frac):
        frac = np.array(frac.cpu())
        t = np.log((1 - frac) / frac)  # Threshold
        ppm[ppm < 1e-12] = 1e-12  # handles small probabilities
        spec = np.log(ppm) - np.log(ppm_bg)  # spec matrix
        spec_revcomp = spec[::-1, ::-1]
        L, W = ppm.shape
        for i in range(0, len(seqs_onehot), 1):
            s = seqs_onehot[i]  # grab the one hot coded sequence
            seqlen = s.shape[1]
            if seqlen < W:  # leave short sequences alone
                continue
            indices = np.arange(seqlen - W + 1)
            conv_signal = signal.correlate2d(spec, s, 'valid')[0]
            seq_motif_sites = indices[conv_signal > t]
            if self.revcomp:
                conv_signal_revcomp = signal.correlate2d(spec_revcomp, s, 'valid')[0]
                seq_motif_sites_revcomp = indices[conv_signal_revcomp > t]
                seq_motif_sites = np.concatenate((seq_motif_sites, seq_motif_sites_revcomp))
            for motif_site in seq_motif_sites:
                s[:, motif_site:motif_site+W] = 0
        seqs = sequences.decode(seqs_onehot, self.alpha)
        return seqs 

def _erase_seqs_containing_motifs(self, seqs_onehot, ppm, ppm_bg, frac):
        frac = np.array(frac.cpu())
        t = np.log((1 - frac) / frac)  # Threshold
        ppm[ppm < 1e-12] = 1e-12  # handles small probabilities
        spec = np.log(ppm) - np.log(ppm_bg)  # spec matrix
        spec_revcomp = spec[::-1, ::-1]
        L, W = ppm.shape
        seqs_onehot_filtered = []
        for i in range(0, len(seqs_onehot), 1):
            s = seqs_onehot[i]  # grab the one hot coded sequence
            if s.shape[1] < W:  # leave short sequences alone
                seqs_onehot_filtered.append(s)
                continue
            conv_signal = signal.correlate2d(spec, s, 'valid')[0]
            s_has_motif = (conv_signal > t).any()
            if self.revcomp:
                conv_signal_revcomp = signal.correlate2d(spec_revcomp, s, 'valid')[0]
                s_has_motif = s_has_motif or (conv_signal_revcomp > t).any()
            if not s_has_motif:
                seqs_onehot_filtered.append(s)
        seqs = sequences.decode(seqs_onehot_filtered, self.alpha)
        return seqs 

def test_consistency_correlate_funcs(self):
        # Compare np.correlate, signal.correlate, signal.correlate2d
        a = np.arange(5)
        b = np.array([3.2, 1.4, 3])
        for mode in ['full', 'valid', 'same']:
            assert_almost_equal(np.correlate(a, b, mode=mode),
                                signal.correlate(a, b, mode=mode))
            assert_almost_equal(np.squeeze(signal.correlate2d([a], [b],
                                                              mode=mode)),
                                signal.correlate(a, b, mode=mode))

            # See gh-5897
            if mode == 'valid':
                assert_almost_equal(np.correlate(b, a, mode=mode),
                                    signal.correlate(b, a, mode=mode))
                assert_almost_equal(np.squeeze(signal.correlate2d([b], [a],
                                                                  mode=mode)),
                                    signal.correlate(b, a, mode=mode)) 

def get_pssm_scores(encoded_sequences, pssm):
  """
  Convolves pssm and its reverse complement with encoded sequences
  and returns the maximum score at each position of each sequence.

  Parameters
  ----------
  encoded_sequences: 3darray
       (N_sequences, N_letters, sequence_length, 1) array
  pssm: 2darray
      (4, pssm_length) array

  Returns
  -------
  scores: 2darray
      (N_sequences, sequence_length)
  """
  encoded_sequences = encoded_sequences.squeeze(axis=3)
  # initialize fwd and reverse scores to -infinity
  fwd_scores = np.full_like(encoded_sequences, -np.inf, float)
  rc_scores = np.full_like(encoded_sequences, -np.inf, float)
  # cross-correlate separately for each base,
  # for both the PSSM and its reverse complement
  for base_indx in range(encoded_sequences.shape[1]):
    base_pssm = pssm[base_indx][None]
    base_pssm_rc = base_pssm[:, ::-1]
    fwd_scores[:, base_indx, :] = correlate2d(
        encoded_sequences[:, base_indx, :], base_pssm, mode='same')
    rc_scores[:, base_indx, :] = correlate2d(
        encoded_sequences[:, -(base_indx + 1), :], base_pssm_rc, mode='same')
  # sum over the bases
  fwd_scores = fwd_scores.sum(axis=1)
  rc_scores = rc_scores.sum(axis=1)
  # take max of fwd and reverse scores at each position
  scores = np.maximum(fwd_scores, rc_scores)
  return scores 

def get_pssm_scores(encoded_sequences, pssm):
  """
  Convolves pssm and its reverse complement with encoded sequences
  and returns the maximum score at each position of each sequence.

  Parameters
  ----------
  encoded_sequences: 3darray
        (num_examples, 1, 4, seq_length) array
  pssm: 2darray
      (4, pssm_length) array

  Returns
  -------
  scores: 2darray
      (num_examples, seq_length) array
  """
  encoded_sequences = encoded_sequences.squeeze(axis=1)
  # initialize fwd and reverse scores to -infinity
  fwd_scores = np.full_like(encoded_sequences, -np.inf, float)
  rc_scores = np.full_like(encoded_sequences, -np.inf, float)
  # cross-correlate separately for each base,
  # for both the PSSM and its reverse complement
  for base_indx in range(encoded_sequences.shape[1]):
    base_pssm = pssm[base_indx][None]
    base_pssm_rc = base_pssm[:, ::-1]
    fwd_scores[:, base_indx, :] = correlate2d(
        encoded_sequences[:, base_indx, :], base_pssm, mode='same')
    rc_scores[:, base_indx, :] = correlate2d(
        encoded_sequences[:, -(base_indx + 1), :], base_pssm_rc, mode='same')
  # sum over the bases
  fwd_scores = fwd_scores.sum(axis=1)
  rc_scores = rc_scores.sum(axis=1)
  # take max of fwd and reverse scores at each position
  scores = np.maximum(fwd_scores, rc_scores)
  return scores 

def _search_minimum_distance(self, ref, buff):
        if len(ref) < self.fl:
            ref = np.r_[ref, np.zeros(self.fl - len(ref))]

        # slicing and windowing one sample by one
        buffmat = view_as_windows(buff, self.fl) * self.win
        refwin = np.array(ref * self.win).reshape(1, self.fl)
        corr = correlate2d(buffmat, refwin, mode='valid')

        return np.argmax(corr) - self.sl 

def test_consistency_correlate_funcs(self):
        # Compare np.correlate, signal.correlate, signal.correlate2d
        a = np.arange(5)
        b = np.array([3.2, 1.4, 3])
        for mode in ['full', 'valid', 'same']:
            assert_almost_equal(np.correlate(a, b, mode=mode),
                                signal.correlate(a, b, mode=mode))
            assert_almost_equal(np.squeeze(signal.correlate2d([a], [b],
                                                              mode=mode)),
                                signal.correlate(a, b, mode=mode))


# Create three classes, one for each complex data type. The actual class
# name will be TestCorrelateComplex###, where ### is the number of bits. 

def test_invalid_shapes(self):
        # By "invalid," we mean that no one
        # array has dimensions that are all at
        # least as large as the corresponding
        # dimensions of the other array. This
        # setup should throw a ValueError.
        a = np.arange(1, 7).reshape((2, 3))
        b = np.arange(-6, 0).reshape((3, 2))

        assert_raises(ValueError, signal.correlate2d, *(a, b), **{'mode': 'valid'})
        assert_raises(ValueError, signal.correlate2d, *(b, a), **{'mode': 'valid'}) 

def test_complex_input(self):
        assert_equal(signal.correlate2d([[1]], [[2j]]), -2j)
        assert_equal(signal.correlate2d([[2j]], [[3j]]), 6)
        assert_equal(signal.correlate2d([[3j]], [[4]]), 12j) 

def getIns(self):
        """Collape matrix into insertions.  Will reduce span on chromosome"""
        pattern = np.zeros((self.upper-self.lower,self.upper + (self.upper-1)%2))
        mid = self.upper/2
        for i in range(self.lower,self.upper):
            pattern[i-self.lower,mid+(i-1)/2]=1
            pattern[i-self.lower,mid-(i/2)]=1
        ins = signal.correlate2d(self.mat,pattern,mode="valid")[0]
        insertion_track = InsertionTrack(self.chrom,self.start + pattern.shape[1]/2, self.end - (pattern.shape[1]/2))
        insertion_track.assign_track(ins)
        return insertion_track 

def forward(ctx, input, filter):
        input, filter = input.detach(), filter.detach()  # detach so we can cast to NumPy
        result = correlate2d(input.numpy(), filter.detach().numpy(), mode='valid')
        ctx.save_for_backward(input, filter)
        return input.new(result) 

def predict_bounding_boxes(self, pointcloud, set_next_bounding_boxes=False, bounding_factor=.5, eps=0.1):
		next_bounding_boxes = []
		for bounding_box in self.bounding_boxes:
			filtered_pointcloud = pointcloud.filter_points()
			filtered_pointcloud_indices = bounding_box.filter_points(filtered_pointcloud, 
																	 bounding_factor=bounding_factor)
			filtered_points = filtered_pointcloud.points[filtered_pointcloud_indices,:]
			x, y = filtered_points[:,0], filtered_points[:,1]
			z = filtered_pointcloud.intensities[filtered_pointcloud_indices]

			# fig = plt.figure()
			# ax = fig.add_subplot(111, projection='3d')
			# ax.scatter(x, y, z)
			# plt.show()

			sorted_x, sorted_y = np.sort(x), np.sort(y)

			resolution = max(eps, min(np.min(np.ediff1d(sorted_x)), np.min(np.ediff1d(sorted_y))))
			h, w = int((np.max(x) - np.min(x)) // resolution) + 1, int((np.max(y) - np.min(y)) // resolution) + 1
			print(h, w, resolution)

			im = -np.ones((h, w)) * 5e-2
			quantized_x = ((filtered_points[:,0] - np.min(x)) // resolution).astype(int)
			quantized_y = ((filtered_points[:,1] - np.min(y)) // resolution).astype(int)
			im[quantized_x, quantized_y] = 1

			mask_h = int(bounding_box.width // resolution) + 1 
			mask_w = int(bounding_box.length // resolution) + 1

			mask = np.ones((mask_h, mask_w))


			# plt.scatter(x, y)
			# plt.show()


			print("mask shape: ", mask.shape)
			cc = signal.correlate2d(im, mask, mode="same")
			center = (np.array([np.argmax(cc) // w, np.argmax(cc) % w]) * resolution +
						np.array([np.min(x), np.min(y)]))
			upper_right = center + np.array([bounding_box.width / 2, bounding_box.length / 2])
			lower_left = center - np.array([bounding_box.width / 2, bounding_box.length / 2])
			theta = bounding_box.angle
			box_pointcloud = PointCloud(np.vstack((upper_right, lower_left)))
			corners = box_pointcloud.rigid_transform(theta, center) + center
			next_bounding_boxes.append([corners.tolist(), theta])
			print(np.argmax(cc) // w, np.argmax(cc) % w, np.argmax(cc), np.max(cc), cc[np.argmax(cc) // w, np.argmax(cc) % w])
			# plt.subplot(1,2,1)
			# plt.imshow(im, cmap='Greys',  interpolation='nearest')
			# plt.subplot(1,2,2)
			# plt.imshow(cc, cmap='Greys',  interpolation='nearest')
			
			
			# plt.show()
		return next_bounding_boxes