Incorporating skycleaver functionality to skyweaver

.gitignore

@@ -23,4 +23,6 @@ eggs/
 core*
 report*
 CMakeCache.txt
+.nextflow
+*.bdm.db
 

Dockerfile.skyweavercpp

@@ -60,7 +60,7 @@ WORKDIR /usr/src/skyweaver
 COPY . . 
 RUN cmake -S . -B build/ -DARCH=native -DPSRDADA_INCLUDE_DIR=/usr/local/include/psrdada \ 
     -DPSRDADACPP_INCLUDE_DIR=/usr/local/include/psrdada_cpp -DSKYWEAVER_NANTENNAS=64 \
-    -DSKYWEAVER_NBEAMS=128 -DSKYWEAVER_NCHANS=64 -DSKYWEAVER_IB_SUBTRACTION=1 -DBUILD_SUBMODULES=1 \
-    -DENABLE_TESTING=1 -DENABLE_BENCHMARK=1 &&\
+    -DSKYWEAVER_NBEAMS=800 -DSKYWEAVER_NCHANS=64 -DSKYWEAVER_IB_SUBTRACTION=1 -DBUILD_SUBMODULES=1 \
+    -DENABLE_TESTING=1 -DENABLE_BENCHMARK=1 -DSKYWEAVER_CB_TSCRUNCH=4 -DSKYWEAVER_IB_TSCRUNCH=4 -DSKYWEAVER_CB_FSCRUNCH=1 -DSKYWEAVER_IB_FSCRUNCH=1 &&\
     make -C build/ -j 16 && make -C build/ install
 

Dockerfile.skyweaverpy

@@ -20,7 +20,8 @@ RUN apt-get update && apt-get install -y \
 
 # Install Python dependencies
 RUN pip install --upgrade pip && \
-    pip install pytest
+    pip install pytest && \
+    pip install pandas
 
 
 # Copy the rest of your application code into the container

README.md

@@ -1,2 +1,177 @@
 # skyweaver
 Implementation of an offline FBFUSE beamformer for the MeerKAT telescope
+
+# Installation
+
+It is easiest to use the software inside a Docker container.  Two dockerfiles is included part of this repository: To compile the c++ and python parts respectively. 
+
+# Usage
+
+## Step 1:  Get delays for the beamformer
+
+### Start with the skyweaver python CLI 
+
+```bash
+alias sw="python /path/to/python/skyweaver/cli.py"
+sw --help
+```
+This will print:
+```console
+usage: skyweaver [-h] {metadata,delays} ...
+
+positional arguments:
+  {metadata,delays}  sub-command help
+    metadata         Tools for observation metadata files
+    delays           Tools for delay files
+
+optional arguments:
+  -h, --help         show this help message and exit
+```
+
+### Get the metadata for the corresponding observation 
+
+This is done outside this repository
+
+### Obtain metadata information
+
+```bash
+sw metadata show <metadata_file>
+```
+
+This will produce an output like the following:
+```console
+sw metadata show bvrmetadata_2024-02-16T10\:50\:46_72275.hdf5
+---------------Array configuration----------------
+Nantennas:          57
+Subarray:           m000,m002,m003,m004,m005,m007,m008,m009,m010,m011,m012,m014,m015,m016,
+                    m017,m018,m019,m020,m021,m022,m023,m024,m025,m026,m027,m029,m030,m031,
+                    m032,m033,m034,m035,m036,m037,m038,m039,m040,m041,m042,m043,m044,m045,
+                    m046,m048,m049,m050,m051,m053,m054,m056,m057,m058,m059,m060,m061,m062,
+                    m063
+Centre frequency:   1284000000.0 Hz
+Bandwidth:          856000000.0 Hz
+Nchannels:          4096
+Sync epoch (UNIX):  1708039531.0
+Project ID:         -
+Schedule block ID:  -
+CBF version:        cbf_dev
+--------------------Pointings---------------------
+#0   J1644-4559          2024-02-16T11:16:08.957000000 until 2024-02-16T11:21:04.865000000 (UTC)
+                         1708082168.957 until 1708082464.865 (UNIX)
+                         72996182384029 until 73502776880016 (SAMPLE CLOCK)
+#1   J1644-4559_Offset1  2024-02-16T11:21:22.092000000 until 2024-02-16T11:26:15.682000000 (UTC)
+                         1708082482.092 until 1708082775.682 (UNIX)
+                         73532269504013 until 74034895583866 (SAMPLE CLOCK)
+#2   J1644-4559_Offset2  2024-02-16T11:26:34.536000000 until 2024-02-16T11:31:25.723000000 (UTC)
+                         1708082794.536 until 1708083085.723 (UNIX)
+                         74067173632022 until 74565685776084 (SAMPLE CLOCK)
+#3   M28                 2024-02-16T11:31:52.503000000 until 2024-02-16T12:01:51.651000000 (UTC)
+                         1708083112.503 until 1708084911.651 (UNIX)
+                         74611533136035 until 77691674512039 (SAMPLE CLOCK)
+#4   J0437-4715          2024-02-16T12:03:06.117000000 until 2024-02-16T12:08:04.364000000 (UTC)
+                         1708084986.117 until 1708085284.364 (UNIX)
+                         77819160304176 until 78329759168140 (SAMPLE CLOCK)
+```
+
+### Create a config file in .yml format
+Here is an example - there are comments to explain each parameter. 
+
+```.yml 
+created_by: Vivek
+beamformer_config:
+  # The total number of beams to be produced (must be a multiple of 32). This needs to be <= the number that SKYWEAVER is compiled for.
+  total_nbeams: 800
+  # The number of time samples that will be accumulated after detection, inside the beamformer
+  tscrunch: 4
+  # The number of frequency channels that will be accumulated after detection, inside the beamformer
+  # Will be coerced to 1 if coherent dedispersion is specified.
+  fscrunch: 1
+  # The Stokes product to be calculated in the beamformer (I=0, Q=1, U=2, V=3)
+  stokes_mode: 0
+  # Enable CB-IB subtraction in the beamformer
+  subtract_ib: True
+
+  # Dispersion measure for coherent / incoherent dedispersion in pc cm^-3
+  # A dispersion plan definition string "
+  #           "(<coherent_dm>:<start_incoherent_dm>:<end_incoherent_dm>:<dm_step>:<tscrunch>) or "
+  #           "(<coherent_dm>:<tscrunch>) "
+  #           "or (<coherent_dm>)")
+# Each DD plan is a "Stream" with zero indexed stream-ids 
+
+ddplan:
+  - "478.6:478.6:478.6:1:1" #stream-id=0
+  - "0.00:478.6:478.6:1:1" #stream-id=1
+
+# every beamset can contain arbitrary set of antennas, corresponding targeted beams, and tiled beams
+# total number of beams across all beamsets should be <= the number of beams that SKYWEAVER is compiled for.
+beam_sets:
+
+  - antenna_set: ['m000','m002','m003','m004','m005','m007','m008','m009','m010','m011',
+                  'm012','m014','m015','m016','m017','m018','m019','m020','m021','m022',
+                  'm023','m024','m025','m026','m027','m029','m030','m031','m032','m033',
+                  'm034','m035','m036','m037','m038','m039','m040','m041','m042','m043',
+                  'm044','m045','m046','m048','m049','m050','m051','m053','m054','m056',
+                  'm057','m058','m059','m060','m061','m062','m063']
+    beams: []
+    tilings:
+      - nbeams: 32
+        reference_frequency: null
+        target: "J1644-4559,radec,16:44:49.273,-45:59:09.71"
+        overlap: 0.9
+```
+
+
+
+### Create delay file for the corresponding pointing
+
+```bash
+sw delays create --pointing-idx 0 --outfile J1644-4559_pointing_0.delays --step 4 bvrmetadata_2024-02-16T10\:50\:46_72275.hdf5 J1644-4559_boresight.yaml
+```
+
+This produces a `.delays` file used for beamforming, and a `.targets` file that contains beam metadata. There are also other files produced here for reproducibility and for visualisation. 
+
+## Step 2: Initialise input and compile skyweaver
+
+### Create a list of dada files that correspond to the pointing
+
+```console
+ls /b/u/vivek/00_DADA_FILES/J1644-4559/2024-02-16-11\:16\:08/L/48/*dada -1 > /bscratch/vivek/skyweaver_tests/J1644-4559_boresight_dadafiles.list
+```
+
+### Compile skyweaver
+
+This is done inside the dockerfile too. Either edit that to produce a docker image that has the software precompiled, or compile separately. 
+
+```bash
+    cmake -S . -B $cmake_tmp_dir -DENABLE_TESTING=0 -DCMAKE_INSTALL_PREFIX=$install_dir -DARCH=native -DPSRDADA_INCLUDE_DIR=/usr/local/include/psrdada -DPSRDADACPP_INCLUDE_DIR=/usr/local/include/psrdada_cpp -DSKYWEAVER_NANTENNAS=64 -DSKYWEAVER_NBEAMS=${nbeams} -DSKYWEAVER_NCHANS=64 -DSKYWEAVER_IB_SUBTRACTION=1 -DCMAKE_BUILD_TYPE=RELEASE -DSKYWEAVER_CB_TSCRUNCH=${tscrunch} -DSKYWEAVER_IB_TSCRUNCH=${tscrunch};
+    cd $cmake_tmp_dir
+    make -j 16
+```
+This compilation produces two binaries: `skyweavercpp` and `skycleaver`
+## Step 3: Run the beamformer
+
+```bash
+
+/path/to/skyweavercpp --input-file J1644-4559_boresight_dadafiles.list --delay-file J1644-4559_pointing_0.delays --output-dir=/bscratch/vivek/skyweaver_out --gulp-size=32768 --log-level=warning --output-level=12 --stokes-mode I
+```
+
+Change the output level to 7 for bright pulsars like J1644-4559. 
+
+This will produce `.tdb` files for the corresponding bridge. Run Step 3 for ALL 64 bridges with their corresponding dada file lists. These are DADA format files with the dimensions of TIME, INCOHERENT DM and BEAM as the order. For stokes I mode, The datatype is `int8_t`. For IQUV it is `char4`. 
+
+## Steo 4: Cleave all bridges to form Filterbanks
+
+Here we cleave the 64 TDB[I/Q/U/V/IV/QU/IQUV] files to produce `NDM*NBEAMS*NSTOKES` number of T(F=64) files. 
+
+to run this, do
+
+```bash
+
+/path/to/skycleaver -r /bscratch/vivek/skyweaver_out --output-dir /bscratch/vivek/skycleaver_out --nsamples-per-block 65536 --nthreads 32 --stream-id 0 --targets_file/bscratch/vivek/skyweaver_out/swdlays_J1644-4559.targets --out-stokes I --required_beams 0 
+
+```
+
+This will produce a standard sigproc format `.fil` file that can be used for traditional processing. 
+
+
+

cmake/compiler_settings.cmake

@@ -5,6 +5,8 @@ if (NOT CMAKE_BUILD_TYPE)
     set(CMAKE_BUILD_TYPE "RELEASE")
 endif ()
 
+#set(CMAKE_VERBOSE_MAKEFILE 1)
+set (CMAKE_CXX_STANDARD 20)
 # Set compiler flags
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC -fopenmp")
 

cpp/skyweaver/CMakeLists.txt

@@ -18,6 +18,8 @@ set(skyweaver_src
    src/Timer.cpp
    src/Transposer.cu
    src/WeightsManager.cu
+   src/SigprocHeader.cpp
+
    )
 
 set(skyweaver_inc
@@ -38,6 +40,10 @@ set(skyweaver_inc
     Timer.hpp
     Transposer.cuh
     WeightsManager.cuh
+    MultiFileWriter.cuh
+    SigprocHeader.hpp
+    SkyCleaver.hpp
+    SkyCleaverConfig.hpp
    )
 
 set(SKYWEAVER_LIBRARIES ${CMAKE_PROJECT_NAME} ${DEPENDENCY_LIBRARIES})
@@ -106,6 +112,13 @@ target_link_libraries(skyweavercpp
     OpenMP::OpenMP_CXX)
 install(TARGETS skyweavercpp DESTINATION bin)
 
+add_executable(skycleaver src/skycleaver_cli.cpp)
+target_link_libraries(skycleaver 
+    ${SKYWEAVER_LIBRARIES} 
+    ${DEPENDENCY_LIBRARIES} 
+    OpenMP::OpenMP_CXX)
+install(TARGETS skycleaver DESTINATION bin)
+
 #install (TARGETS ... DESTINATION bin)
 #install (TARGETS ${CMAKE_PROJECT_NAME}
 #         RUNTIME DESTINATION bin
@@ -126,4 +139,4 @@ if(ENABLE_BENCHMARK)
     COMMAND beamformer_bench --benchmark_counters_tabular=true
     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
   add_subdirectory(benchmark)
-endif()
+endif()

cpp/skyweaver/DescribedVector.hpp

@@ -127,6 +127,23 @@ struct DescribedVector {
         _vector.resize(calculate_nelements());
     }
 
+    /**
+     * @brief Construct a new Described Vector object of specific size
+     *
+     * @param sizes The sizes of the dimensions (must match the number of
+     * dimensions)
+     */
+    DescribedVector(std::initializer_list<std::size_t> sizes, value_type default_value)
+        : _dms_stale(true), _frequencies_stale(true), _sizes(sizes),
+          _dims{dims...}, _tsamp(0.0)
+    {
+        if(_sizes.size() != sizeof...(dims)) {
+            throw std::invalid_argument(
+                "Number of sizes must match number of dimensions");
+        }
+        _vector.resize(calculate_nelements(), default_value);
+    }
+
     /**
      * @brief Destroy the Described Vector object
      *
@@ -213,6 +230,12 @@ struct DescribedVector {
      */
     auto const& operator[](std::size_t idx) const { return _vector[idx]; }
 
+
+    // also the at() method
+    auto& at(std::size_t idx) { return _vector[idx]; }
+    auto const& at(std::size_t idx) const { return _vector[idx]; }
+
+
     /**
      * @brief Resize the dimensions of the vector
      *
@@ -343,12 +366,13 @@ struct DescribedVector {
      * @param freqs
      */
     void frequencies(FrequenciesType const& freqs)
-    {
+    {        
         if(freqs.size() != get_dim_extent<FreqDim>()) {
             throw std::runtime_error("Invalid number of frequecies passed.");
         }
         _frequencies_stale = false;
         _frequencies       = freqs;
+
     }
 
     /**
@@ -363,6 +387,7 @@ struct DescribedVector {
         }
         _frequencies_stale = false;
         _frequencies.resize(1, freq);
+        _frequencies[0] = freq;
     }
 
     /**
@@ -648,7 +673,7 @@ using BTFPowersH = DescribedVector<thrust::host_vector<T, PinnedAllocator<T>>,
 template <typename T>
 using BTFPowersD =
     DescribedVector<thrust::device_vector<T>, BeamDim, TimeDim, FreqDim>;
-// Incoherent dedisperser outputs
+// Incoherent dedisperser outputs, also used for skycleaver
 template <typename T>
 using TDBPowersH = DescribedVector<thrust::host_vector<T, PinnedAllocator<T>>,
                                    TimeDim,
@@ -667,6 +692,19 @@ template <typename T>
 using FPAStatsD =
     DescribedVector<thrust::device_vector<T>, FreqDim, PolnDim, AntennaDim>;
 
+
+//skycleaver vectors
+template <typename T>
+using TDBPowersStdH = DescribedVector<std::vector<T>,
+                                   TimeDim,
+                                   DispersionDim,
+                                   BeamDim,
+                                   PolnDim>;
+template <typename T>
+using TFPowersStdH = DescribedVector<std::vector<T>,
+                                  TimeDim,
+                                  FreqDim>;
+
 } // namespace skyweaver
 
 #endif // SKYWEAVER_DESCRIBEDVECTORS_HPP

cpp/skyweaver/FileOutputStream.hpp

@@ -48,6 +48,8 @@ class FileOutputStream
         std::string _full_path;
         std::size_t _bytes_requested;
         std::size_t _bytes_written;
+        std::string _temporary_suffix;
+        std::string _temporary_path;
         std::ofstream _stream;
     };
 

cpp/skyweaver/Header.hpp

@@ -37,6 +37,26 @@ class Header
     template <typename T>
     T get(char const* key) const;
 
+
+    /**
+      * @brief      Get a value from the header or return a default value
+      *
+      * @param      key            An ASCII key (the name of the value to get)
+      * @param[in]  default_value  The default value to return if the key is not found
+      *
+      * @tparam     T              The data type of the parameter to be read
+      *
+      * @return     The value corresponding to the given key or the default value
+     */
+    template <typename T>
+    T get_or_default(char const* key, T default_value) const{
+        if(has_key(key)){
+            return get<T>(key);
+        }
+        BOOST_LOG_TRIVIAL(warning) << "Header did not have key: " << key << ". Returning default value: " << default_value;
+        return default_value;
+    }
+
     /**
      * @brief      Set a value in the header
      *
@@ -72,6 +92,7 @@ class Header
 
   private:
     void fetch_header_string(char const* key) const;
+    bool has_key(char const* key) const;
 
   private:
     psrdada_cpp::RawBytes& _header;