Small speed up

Add normalised value

gsrp_smart_util.py

@@ -90,9 +90,11 @@ def mul(mem1, mem2, cc, t_max, id1, id2, n_ind):
     # assume len(id2) == 1
     idx1, idx2 = np.meshgrid(np.arange(len(mem1[0])), np.arange(len(mem2[0])))
     idx1, idx2 = idx1.ravel(), idx2.ravel()
-    out_tij = np.concatenate((mem1[1][:, idx1], mem2[1][:, idx2]), axis=0)
-    mask = (np.abs(out_tij[:-1] - out_tij[-1:]) <= t_max).all(0)
-    out_tij = out_tij[:, mask]
+    out_tij = np.empty((len(id1) + 1, len(idx1)), mem1[1].dtype)
+    np.take(mem1[1], idx1, axis=1, out=out_tij[:-1])
+    np.take(mem2[1], idx2, axis=1, out=out_tij[-1:])
+    mask = reduce_all(np.abs(out_tij[:-1] - out_tij[-1:]) <= t_max) # Faster than numpy for large array
+    out_tij = np.compress(mask, out_tij, axis=1)
     idx1, idx2 = idx1[mask], idx2[mask]
     out_val = mem1[0][idx1] * mem2[0][idx2]
     tij_dep = out_tij[:-1] - out_tij[-1:]
@@ -104,12 +106,19 @@ def mul(mem1, mem2, cc, t_max, id1, id2, n_ind):
 
 def mask_val(mem, val):
     mask = mem[0] >= val
-    return mem[0][mask], mem[1][:, mask]
+    return mem[0][mask], np.compress(mask, mem[1], axis=1)
+
+
+def reduce_all(arr):
+    out = arr[0].copy()
+    for a in arr[1:]:
+        out &= a
+    return out
 
 
 def mask_lim(mem, tij_min, tij_max, t_max):
-    mask = ((tij_max - t_max <= mem[1]) & (mem[1] <= tij_min + t_max)).all(0)
-    return mem[0][mask], mem[1][:, mask]
+    mask = reduce_all((tij_max - t_max <= mem[1]) & (mem[1] <= tij_min + t_max))  # Faster than numpy for large array
+    return mem[0][mask], np.compress(mask, mem[1], axis=1)
 
 
 def constrained_argmax(mem, cc, tij_ind, curr_tij, used_tij, t_max, n_ind):
@@ -134,7 +143,7 @@ def smart_gsrp(cc, n_ind, n_tot, t_max, tree, program, clean_list):
             if op == 'mem':
                 if i == 0:
                     memory[(0, j)] = mask_val((np.concatenate((cc[op.right, -t_max:], cc[op.right, :t_max + 1])),
-                                               np.arange(-t_max, t_max + 1)[np.newaxis]), val)
+                                               np.arange(-t_max, t_max + 1, dtype=np.int32)[np.newaxis]), val)
                 else:
                     memory[(i, j)] = mask_val(memory[(op.left, op.right)], val)
             else:  # op == 'mul'
@@ -161,10 +170,4 @@ def smart_gsrp(cc, n_ind, n_tot, t_max, tree, program, clean_list):
         for p in clean_list[i]:
             del memory[p]
 
-    return val, tij
-
-
-
-
-
-
+    return np.log10(val), tij

gsrp_tdoa_hyperres.py

@@ -78,10 +78,10 @@ def corr(data, pos, w_size, max_tdoa, decimate=1, mode='prepare', hyper=True):
     else:
         raise ValueError(f'Unknown mode {mode}')
 
-    tdoas = np.zeros((len(pos), num_channel_pairs + 1), np.float32)
+    tdoas = np.zeros((len(pos), num_channel_pairs + 2), np.float32)
 
     if hyper:   # prepare hyper res
-        tdoas2 = np.zeros((len(pos), num_channel_pairs + 1), np.float32)
+        tdoas2 = np.zeros((len(pos), num_channel_pairs + 2), np.float32)
         poly = PolynomialFeatures(2)
         lin = LinearRegression()
         pipe = Pipeline([('poly', poly), ('lin', lin)])
@@ -98,8 +98,8 @@ def corr(data, pos, w_size, max_tdoa, decimate=1, mode='prepare', hyper=True):
                             ).named_steps['lin'].coef_
             der = np.zeros((num_channels - 1, num_channels - 1))
             der[ind] = coef[4:]
-            poly_min = np.linalg.lstsq(der + der.T, -coef[1:4], rcond=None)[0] + mean
-            return pipe.predict(poly_min[np.newaxis]).item(), mat @ poly_min
+            poly_min = np.linalg.lstsq(der + der.T, -coef[1:4], rcond=None)[0]
+            return np.log10(pipe.predict(poly_min[np.newaxis]).item()), mat @ (poly_min + mean)
 
     cc = np.empty((num_channel_pairs, dw_size), np.float32)
     for i in trange(len(pos)):
@@ -109,25 +109,27 @@ def corr(data, pos, w_size, max_tdoa, decimate=1, mode='prepare', hyper=True):
         fft = np.asarray(fft, dtype=np.complex64)
         cc[:] = irfft(fft[:, v2] * np.conj(fft[:, v1]), axis=0).T
         cc -= cc.min(-1, keepdims=True)
+        maxs = cc.max(1, keepdims=True)
+        cc /= maxs
+        maxs = np.log10(maxs.prod())
         if mode == 'prepare':
-            tdoas[i, 0], index = c_corr.c_corr_at(cc, tausf)
-            tdoas[i, 1:] = ((tausf[:, index] + dw_size // 2) % dw_size) - dw_size // 2
+            tdoas[i, :2], index = c_corr.c_corr_at(cc, tausf)
+            tdoas[i, 2:] = ((tausf[:, index] + dw_size // 2) % dw_size) - dw_size // 2
         elif mode == 'on-the-fly':
-            tdoas[i, 0], tdoas[i, 1:] = c_corr.c_corr_all(cc, cc_size//2, num_channels - 1)
+            tdoas[i, :2], tdoas[i, 2:] = c_corr.c_corr_all(cc, cc_size//2, num_channels - 1)
         elif mode == 'smart':
-            maxs = cc.max(1, keepdims=True)
-            cc /= maxs
-            val, tdoas[i, 1:] = smart_gsrp(cc, num_channels - 1, num_channel_pairs, cc_size // 2,
+            tdoas[i, :2], tdoas[i, 2:] = smart_gsrp(cc, num_channels - 1, num_channel_pairs, cc_size // 2,
                                            tree, program, clean_list)
-            tdoas[i, 0] = np.log10(val * maxs.prod())
-            cc *= maxs
         else:
             raise ValueError(f'Unknown mode {mode}')
+        tdoas[i, 1] += maxs
 
         if hyper:
-            tdoas2[i, 0], tdoas2[i, 1:] = _hyperres(tdoas[i, 1:], cc)
-
+            tdoas2[i, :2], tdoas2[i, 2:] = _hyperres(tdoas[i, 2:], cc)
+            tdoas2[i, 1] += maxs
+    tdoas[:, :2] *= 20
     if hyper:
+        tdoas2[:, :2] *= 20
         return np.hstack((np.expand_dims(pos, -1), tdoas)), np.hstack((np.expand_dims(pos, -1), tdoas2))
     else:
         return np.hstack((np.expand_dims(pos, -1), tdoas))
@@ -219,8 +221,9 @@ if __name__ == "__main__":
 
     parser.add_argument('infile', type=str, help='The sound file to process.')
     parser.add_argument('outfile', type=str, help='The text or npy file to write results to. Each row gives the '
-                                                  'position (in samples), cross-correlation product, the independent '
-                                                  'TDOAs (in samples), and TDOAs derived from the independent ones. ')
+                                                  'position (in samples), cross-correlation product in decibel '
+                                                  '(normalized and unormalized), the independent TDOAs (in samples), '
+                                                  'and TDOAs derived from the independent ones.')
 
     group1 = parser.add_argument_group('Channels')
     group1.add_argument('-c', '--channels', type=intlist, default=None,