search

UCBerkeleySETI / hyperseti

A SETI / technosignature search code to find intelligent life beyond Earth
https://hyperseti.readthedocs.io
10 stars 4 forks source link

The example in README.md is broken #48

Closed texadactyl closed 2 years ago

texadactyl commented 2 years ago

The example needs a re-vamp. FWIW, the following re-vamp works although the output is probably not 100% desirable.

E.g. the find_et call on a Voyager 1 .h5 file takes 99 seconds as opposed to turbo_seti taking less than 10 seconds. And, the returned dataframe suggests that there are other issues.

# Change the following line to make a file search work.
FILEPATH= "/home/texadactyl/hyperseti/test/test_data/Voyager1.single_coarse.fine_res.h5"

import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
from hyperseti import dedoppler, run_pipeline, find_et
from hyperseti.plotting import imshow_dedopp, imshow_waterfall

# Create a drifting test signal
N_timestep, N_chans = 32, 256
test_data = np.ones(shape=(N_timestep, 1, N_chans ), dtype='float32')
for ii in range(N_timestep):
    test_data[ii, :, N_chans // 2 + ii] = 100.0

# Create basic metadata
metadata = {'frequency_start': 1000*u.MHz, 'time_step': 1.0*u.s, 'frequency_step': 1.0*u.Hz}

# Run dedoppler
dedopp, metadata = dedoppler(test_data,
                             metadata,
                             boxcar_size=1,
                             max_dd=4.0)

# Imshow output
plt.figure(figsize=(8, 3))
plt.subplot(1,2,1)
imshow_waterfall(np.log(test_data), metadata)
plt.subplot(1,2,2)
imshow_dedopp(np.log(dedopp), metadata)
plt.tight_layout()
plt.savefig("./figures.1.png")
print("metadata:", metadata)

# ===================================================
# Can also search for hits in the dedoppler spectra
# ===================================================
from hyperseti import hitsearch
hits = hitsearch(dedopp.data, metadata, threshold=500)

from hyperseti.plotting import overlay_hits
overlay_hits(hits)
plt.savefig("./figures.2.png")

# ===========================================================================
# The dedoppler and hitsearch can be done in one line with run_pipeline(). 
# Data can be boxcar averaged to look for wider-band signals, 
# and to retrieve signal-to-noise for signals with large drift rates
# ===========================================================================

print("run_pipeline .....")
peaks = run_pipeline(test_data, metadata, max_dd=1.0, threshold=100, 
                     n_boxcar=5, merge_boxcar_trials=True)
print("Returned peaks:", peaks)

# ===================================================
# Can also search for hits in the dedoppler spectra
# ===================================================
print("find_et from file {} .....".format(FILEPATH))
dframe = find_et(FILEPATH, filename_out='./hits.csv', gulp_size=2**18, max_dd=4.0, threshold=50)
print("Returned dataframe:", dframe)
texadactyl commented 2 years ago

stdout:

metadata: {'frequency_start': <Quantity 1000. MHz>, 'time_step': <Quantity 1. s>, 'frequency_step': <Quantity 1. Hz>, 'drift_rate_start': <Quantity -4. Hz / s>, 'drift_rate_step': <Quantity 0.03125 Hz / s>, 'drift_rates': array([-4.     , -3.96875, -3.9375 , -3.90625, -3.875  , -3.84375,
       -3.8125 , -3.78125, -3.75   , -3.71875, -3.6875 , -3.65625,
       -3.625  , -3.59375, -3.5625 , -3.53125, -3.5    , -3.46875,
       -3.4375 , -3.40625, -3.375  , -3.34375, -3.3125 , -3.28125,
       -3.25   , -3.21875, -3.1875 , -3.15625, -3.125  , -3.09375,
       -3.0625 , -3.03125, -3.     , -2.96875, -2.9375 , -2.90625,
       -2.875  , -2.84375, -2.8125 , -2.78125, -2.75   , -2.71875,
       -2.6875 , -2.65625, -2.625  , -2.59375, -2.5625 , -2.53125,
       -2.5    , -2.46875, -2.4375 , -2.40625, -2.375  , -2.34375,
       -2.3125 , -2.28125, -2.25   , -2.21875, -2.1875 , -2.15625,
       -2.125  , -2.09375, -2.0625 , -2.03125, -2.     , -1.96875,
       -1.9375 , -1.90625, -1.875  , -1.84375, -1.8125 , -1.78125,
       -1.75   , -1.71875, -1.6875 , -1.65625, -1.625  , -1.59375,
       -1.5625 , -1.53125, -1.5    , -1.46875, -1.4375 , -1.40625,
       -1.375  , -1.34375, -1.3125 , -1.28125, -1.25   , -1.21875,
       -1.1875 , -1.15625, -1.125  , -1.09375, -1.0625 , -1.03125,
       -1.     , -0.96875, -0.9375 , -0.90625, -0.875  , -0.84375,
       -0.8125 , -0.78125, -0.75   , -0.71875, -0.6875 , -0.65625,
       -0.625  , -0.59375, -0.5625 , -0.53125, -0.5    , -0.46875,
       -0.4375 , -0.40625, -0.375  , -0.34375, -0.3125 , -0.28125,
       -0.25   , -0.21875, -0.1875 , -0.15625, -0.125  , -0.09375,
       -0.0625 , -0.03125,  0.     ,  0.03125,  0.0625 ,  0.09375,
        0.125  ,  0.15625,  0.1875 ,  0.21875,  0.25   ,  0.28125,
        0.3125 ,  0.34375,  0.375  ,  0.40625,  0.4375 ,  0.46875,
        0.5    ,  0.53125,  0.5625 ,  0.59375,  0.625  ,  0.65625,
        0.6875 ,  0.71875,  0.75   ,  0.78125,  0.8125 ,  0.84375,
        0.875  ,  0.90625,  0.9375 ,  0.96875,  1.     ,  1.03125,
        1.0625 ,  1.09375,  1.125  ,  1.15625,  1.1875 ,  1.21875,
        1.25   ,  1.28125,  1.3125 ,  1.34375,  1.375  ,  1.40625,
        1.4375 ,  1.46875,  1.5    ,  1.53125,  1.5625 ,  1.59375,
        1.625  ,  1.65625,  1.6875 ,  1.71875,  1.75   ,  1.78125,
        1.8125 ,  1.84375,  1.875  ,  1.90625,  1.9375 ,  1.96875,
        2.     ,  2.03125,  2.0625 ,  2.09375,  2.125  ,  2.15625,
        2.1875 ,  2.21875,  2.25   ,  2.28125,  2.3125 ,  2.34375,
        2.375  ,  2.40625,  2.4375 ,  2.46875,  2.5    ,  2.53125,
        2.5625 ,  2.59375,  2.625  ,  2.65625,  2.6875 ,  2.71875,
        2.75   ,  2.78125,  2.8125 ,  2.84375,  2.875  ,  2.90625,
        2.9375 ,  2.96875,  3.     ,  3.03125,  3.0625 ,  3.09375,
        3.125  ,  3.15625,  3.1875 ,  3.21875,  3.25   ,  3.28125,
        3.3125 ,  3.34375,  3.375  ,  3.40625,  3.4375 ,  3.46875,
        3.5    ,  3.53125,  3.5625 ,  3.59375,  3.625  ,  3.65625,
        3.6875 ,  3.71875,  3.75   ,  3.78125,  3.8125 ,  3.84375,
        3.875  ,  3.90625,  3.9375 ,  3.96875,  4.     ]), 'boxcar_size': 1, 'n_integration': 32, 'integration_time': <Quantity 1. s>, 'obs_len': <Quantity 32. s>, 'output_dims': ('drift_rate', 'feed_id', 'frequency')}
run_pipeline .....
Returned peaks: Empty DataFrame
Columns: []
Index: []
find_et from file /home/texadactyl/hyperseti/test/test_data/Voyager1.single_coarse.fine_res.h5 .....
## TOTAL TIME: 98.93s ##

Returned dataframe:     drift_rate      f_start         snr  driftrate_idx  channel_idx  boxcar_size  beam_idx  n_integration
20    0.095665  8419.921874  825.091370          428.0          0.0         32.0       0.0           16.0
2    -0.200896  8419.189457  126.649879          397.0     262142.0          1.0       0.0           16.0
7     1.033181  8419.921664   80.166428          526.0         75.0          2.0       0.0           16.0
8    -1.281910  8419.189695   62.036343          284.0     262057.0          2.0       0.0           16.0
16   -1.501939  8419.189776   61.668865          261.0     262028.0          4.0       0.0           16.0
9     1.310609  8419.921575   61.271889          555.0        107.0          2.0       0.0           16.0
10    1.454106  8419.921552   60.768215          570.0        115.0          2.0       0.0           16.0
11    1.549771  8419.921533   60.275242          580.0        122.0          2.0       0.0           16.0
13    1.635870  8419.921499   51.660404          589.0        134.0          2.0       0.0           16.0
3    -0.660088  8419.189563   51.307178          349.0     262104.0          1.0       0.0           16.0
4     0.707920  8419.921745   51.276970          492.0         46.0          1.0       0.0           16.0