Merge branch 'spol_fix' into 'devel'

Fix detect accumulate in Dualpol Spectrometer

See merge request mpifr-bdg/psrdada_cpp!16

psrdada_cpp/effelsberg/edd/detail/GatedSpectrometer.cu

@@ -371,16 +371,16 @@ void GatedSpectrometer<HandlerType, InputType, OutputType>::process(SinglePolari
             thrust::raw_pointer_cast(inputDataStream->_channelised_voltage_G0.data()),
             thrust::raw_pointer_cast(outputDataStream->G0.data.a().data()),
             _nchans,
-            inputDataStream->_channelised_voltage_G0.size() / _nchans,
-            _naccumulate / _nBlocks,
+            inputDataStream->_channelised_voltage_G0.size(),
+            _naccumulate,
             1, 0., 1, 0);
 
   kernels::detect_and_accumulate<IntegratedPowerType> <<<1024, 1024, 0, _proc_stream>>>(
             thrust::raw_pointer_cast(inputDataStream->_channelised_voltage_G1.data()),
             thrust::raw_pointer_cast(outputDataStream->G1.data.a().data()),
             _nchans,
-            inputDataStream->_channelised_voltage_G1.size() / _nchans,
-            _naccumulate / _nBlocks,
+            inputDataStream->_channelised_voltage_G1.size(),
+            _naccumulate,
             1, 0., 1, 0);
 
     // count saturated samples

psrdada_cpp/effelsberg/edd/src/GatedSpectrometer.cu

@@ -116,8 +116,12 @@ __global__ void update_baselines(float*  __restrict__ baseLineG0,
        NG1 += stats_G1[i];
     }
 
+    if (NG0 > 0)
         baseLineG0[0] = baseLineNG0[0] / NG0;
+
+    if (NG1 > 0)
         baseLineG1[0] = baseLineNG1[0] / NG1;
+
     baseLineNG0[0] = 0;
     baseLineNG1[0] = 0;
 }
@@ -375,8 +379,14 @@ void GatedPowerSpectrumOutput::swap(cudaStream_t &_proc_stream)
 
 void GatedPowerSpectrumOutput::data2Host(cudaStream_t &_d2h_stream)
 {
-    // copy data to host if block is finished
-  CUDA_ERROR_CHECK(cudaStreamSynchronize(_d2h_stream));
+    BOOST_LOG_TRIVIAL(debug) << "Copying data2Host for " << G0._noOfBitSets.size() << " blocks.";
+
+    // memory layout of output block
+    // --------------------------------------------------------------------------------------------
+    // | spectrum ND0 | spectrum ND1 | Nsamples ND0 | Noverflow ND0 |Nsamples ND1 | Noverflow ND1 |
+    // | spectrum ND0 | spectrum ND1 | Nsamples ND0 | Noverflow ND0 |Nsamples ND1 | Noverflow ND1 |
+    // | ...                                                                                      |
+    // --------------------------------------------------------------------------------------------
 
   for (size_t i = 0; i < G0._noOfBitSets.size(); i++)
   {

psrdada_cpp/effelsberg/edd/test/src/DetectorAccumulatorTester.cu

@@ -139,8 +139,6 @@ TEST_P(detect_and_accumulate32bit, no_stride)
             params.naccumulate,
             1, 0., 1, 0);
 
-    cudaDeviceSynchronize();
-
     thrust::host_vector<float> output_host = output;
 
     for (size_t i =0; i < output.size(); i++)
@@ -173,8 +171,6 @@ TEST_P(detect_and_accumulate32bit, w_stride)
             params.naccumulate,
             1, 0., 2, 0);
 
-    cudaDeviceSynchronize();
-
     thrust::host_vector<float> output_host = output;
 
     for (size_t i =0; i < params.nchan * params.nspectra; i++)

psrdada_cpp/effelsberg/edd/test/src/GatedSpectrometerTest.cu

@@ -496,10 +496,72 @@ TEST_P(ExecutionTests, SinglePolOutput)
     // Test buffer consistency with input parameters
     EXPECT_EQ(bufferLayout.getNHeaps(), params.nHeaps);
 
-    psrdada_cpp::NullSink sink;
+    struct output_pointer
+    {
+        float* spectrum;
+        size_t* nbitsset;
+        size_t* nsaturated;
+        output_pointer(psrdada_cpp::RawBytes& block, size_t _nchans, size_t ND, size_t offset)
+        {
+            size_t specsize = _nchans * sizeof(float);
+            size_t memoffset = 2 * offset * (specsize + 2 * sizeof(size_t));
+
+            spectrum = (float*)static_cast<void *>(block.ptr() + ND * specsize + memoffset);
+            nbitsset = (size_t*)static_cast<void *>(block.ptr() + memoffset + 2 * specsize + ND * 2 * sizeof(size_t));
+            nsaturated =(size_t*)static_cast<void *>(block.ptr() + memoffset + 2 * specsize + ND * 2 * sizeof(size_t) + sizeof(size_t));
+        }
+    };
+
+    class TestSink
+    {
+        size_t _counter;
+        size_t _nchans;
+    public:
+        TestSink(size_t nchans) : _nchans(nchans), _counter(0){};
+        ~TestSink(){};
+        void init(psrdada_cpp::RawBytes&){};
+        bool operator()(psrdada_cpp::RawBytes& block){
+            _counter++;
+            output_pointer ND0(block, _nchans, 0, 0);
+            output_pointer ND1(block, _nchans, 1, 0);
+
+            if (_counter == 1)
+            { // First sepctrum all zeros in ND on and ND off
+                for (size_t i=1; i < _nchans; i++)
+                {
+                    EXPECT_NEAR(ND0.spectrum[i], 0, 1E-10) <<  ", i = " << i;
+                    EXPECT_NEAR(ND1.spectrum[i], 0, 1E-10) <<  ", i = " << i;
+                }
+            }
+            else if (_counter == 2)
+            { // Second spectrum all data in ND1 (and flat)
+
+                for (size_t i=1; i < _nchans; i++)
+                {
+                    EXPECT_GT(ND0.spectrum[i], 1) << i;
+                    EXPECT_NEAR(ND0.spectrum[i], ND0.spectrum[_nchans / 2], 1E-1) << ", i = " << i;
+                    EXPECT_NEAR(ND1.spectrum[i], 0, 1E-10) << i;
+                }
+            }
+            else if (_counter == 3)
+            { // Third spectrum all zeros in ND0
+                for (size_t i=1; i < _nchans; i++)
+                {
+                    EXPECT_NEAR(ND0.spectrum[i], 0, 1E-10) << i;
+                    EXPECT_GT(ND1.spectrum[i], 1) << i;
+                    EXPECT_NEAR(ND1.spectrum[i], ND1.spectrum[_nchans / 2], 1E-1) <<  ", i = " << i;
+                }
+
+            }
+
+            return false;
+        };
+    };
+    TestSink sink(params.fft_length / 2 + 1);
+
 
     psrdada_cpp::effelsberg::edd::GatedSpectrometer<
-        psrdada_cpp::NullSink,
+        TestSink,
         psrdada_cpp::effelsberg::edd::SinglePolarizationInput,
         psrdada_cpp::effelsberg::edd::GatedPowerSpectrumOutput>
             spectrometer(bufferLayout,
@@ -509,15 +571,28 @@ TEST_P(ExecutionTests, SinglePolOutput)
                     sink);
 
     char  *raw_buffer = new char[inputBufferSize];
-    memset(raw_buffer, 0, inputBufferSize);
 
     psrdada_cpp::RawBytes buff(raw_buffer, inputBufferSize, inputBufferSize);
-
     EXPECT_NO_THROW(spectrometer.init(buff)) << params.msg;
-    for (int i = 0; i < 5; i++)
-    {
+    // 1 Block, all zeros
+    memset(raw_buffer, 0, inputBufferSize);
+    EXPECT_NO_THROW(spectrometer(buff)) << params.msg;
+
+    // 2 Block, set value
+    raw_buffer[122] =  int8_t(13);
+    EXPECT_NO_THROW(spectrometer(buff)) << params.msg;
+
+    // 3 Block, set value and SCI
+    raw_buffer[122] =  int8_t(13);
+    uint64_t* sc_items = reinterpret_cast<uint64_t*>(raw_buffer + params.nHeaps * heapSize);
+    for (size_t i = 0; i < params.nHeaps; i++)
+        SET_BIT(sc_items[i], 0);
+    EXPECT_NO_THROW(spectrometer(buff)) << params.msg;
+
+    // Additional three executions to get blocks out
+    EXPECT_NO_THROW(spectrometer(buff)) << params.msg;
+    EXPECT_NO_THROW(spectrometer(buff)) << params.msg;
     EXPECT_NO_THROW(spectrometer(buff)) << params.msg;
-    }
 
     delete [] raw_buffer;
 }