Ares.c

/******************************************************************************************************************************************************
 * Ares.c
 *
 * Software to read beam data from (TEL_SHM) and create the FRB_SHM.
 * COMPILE: gcc -g -O3 -Wall -Wextra -std=gnu99 Ares.c -o Ares -lm -lc
 * USAGE: ./Ares [observation duration (seconds) = 3600]
 mpicc -std=c99 -mtune=nocona -O3 -ffast-math -msse2 -lm -I. -g -fopenmp -O4
 -lgomp -I /lustre_archive/apps/correlator/cuda-9.1/include -o Ares Ares.c
 -D_XOPEN_SOURCE=700
 *
 *
 ******************************************************************************************************************************************************/

#define _POSIX_C_SOURCE 200809L
#include "acqpsr.h"
#include "frb_shm.h"
#include "gmrt_newcorr.h"
#include <float.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/shm.h>
#include <sys/time.h>
#include <time.h>

// #define _XOPEN_SOURCE 700
// #define _DEFAULT_SOURCE
#include <unistd.h>

#define BeamHdrKey 1050

int commSize, myrank;
int bandsign;
float bandwidth;

static BeamHeaderType *dataHdr_TEL;
static RecType *Rec;
GlobalInfoType *dataBuffer_TEL;
BeamHeaderType *dataHdr_FRB;
DataBuffer *dataBuffer_FRB;
unsigned char *Shmp, *blkp_TEL, *blkp_FRB;
unsigned int RecNum = 0, TEL_to_FRB_RecNum = 0, RecNum_FRB = 0, DataSeq = 0,
             DataSeq_FRB = 0;
// RecNum goes from 0 to MaxRecs (MaxRecs = 8 for TEL_SHM and 12*32 for
// FRB_SHM.) Only for the sake of this program, I am pretending as if FRB_SHM is
// 12 * 32 blocks long, with each block equal in size to that of TEL_SHM. In the
// FRB_SHM, the curRecord and curBlock are updated after accumulating 32 blocks
// of TEL_SHM. DataSeq keeps track of the number of blocks (or records) read
// since the starting of the observation. It is, in a sense, cummulative RecNum.

unsigned int *shm_marker_p;
double Data_Time[MaxRecs + 1] = {0.0};

int ind = 0, beam_ID = 0;
struct tm *local_t;
char time_string[40];
int usec = 0;
double blk_nano = 0;

struct timeval tv1;
// struct timezone tz1;
typedef struct {
  int tz_minuteswest; /* minutes west of Greenwich */
  int tz_dsttime;     /* type of DST correction */
} timezone1;

timezone1 tz1;
int nBeams;
long int CurrRecSize;
long dataSize;
int curr_chans;
double Freq, chWidth, chan_freq;

// Filter related parameters + initialization
int block_size;
float thresig;
int filltype;
int nchans;
float outmean;
float outstd;
float width;
float tsamp;
size_t raw_size;
size_t nsamples;
size_t block_len;
size_t n_full;
uint8_t *block_ptr;

/*
 * rfi_mitigation_omp_fixed.c
 * C conversion of original C++ code by Raghav Wani (Last edited 18th Feb 2026)
 * Converted to C11 + OpenMP.  Correctness-fixed parallel version.
 *
 * Compile:
 *   gcc -O2 -std=c11 -fopenmp -o rfi_mitigation rfi_mitigation_omp_fixed.c -lm
 *
 * Control thread count at runtime:
 *   OMP_NUM_THREADS=8 ./rfi_mitigation <file> <block_size> <threshold>
 *
 */

/* =========================================================================
 * Portability helpers
 * ====================================================================== */
static float clampf(float v, float lo, float hi) {
  if (v < lo)
    return lo;
  if (v > hi)
    return hi;
  return v;
}
static double now_usec(void) { return omp_get_wtime() * 1e6; }
static double now_msec(void) { return omp_get_wtime() * 1e3; }

/* =========================================================================
 * Error macros
 * ====================================================================== */
#define FATAL(msg)                                                             \
  do {                                                                         \
    fprintf(stderr, "Error: %s\n", msg);                                       \
    exit(1);                                                                   \
  } while (0)
#define FATALF(fmt, ...)                                                       \
  do {                                                                         \
    fprintf(stderr, "Error: " fmt "\n", __VA_ARGS__);                          \
    exit(1);                                                                   \
  } while (0)

/* =========================================================================
 * Kahan compensated summation helper
 * Gives the same result as serial left-to-right double accumulation
 * regardless of the order in which partial sums arrive.
 * ====================================================================== */
typedef struct {
  double sum;
  double c;
} KahanSum;

static inline void kahan_add(KahanSum *ks, double v) {
  double y = v - ks->c;
  double t = ks->sum + y;
  ks->c = (t - ks->sum) - y;
  ks->sum = t;
}

/* =========================================================================
 * Thread-local Box-Muller RNG  (Xorshift64 + Box-Muller)
 * ====================================================================== */
typedef struct {
  unsigned long long state;
} Rng;

static void rng_init(Rng *r, int extra_seed) {
  r->state = (unsigned long long)time(NULL) ^ (unsigned long long)(uintptr_t)r ^
             ((unsigned long long)extra_seed * 6364136223846793005ULL);
  r->state ^= r->state << 13;
  r->state ^= r->state >> 7;
  r->state ^= r->state << 17;
}
static double rng_uniform(Rng *r) {
  r->state ^= r->state << 13;
  r->state ^= r->state >> 7;
  r->state ^= r->state << 17;
  return (double)(r->state >> 11) / (double)(1ULL << 53);
}
static float rng_normal(Rng *r, float mean, float std) {
  double u1 = rng_uniform(r);
  double u2 = rng_uniform(r);
  if (u1 < 1e-300)
    u1 = 1e-300;
  double z = sqrt(-2.0 * log(u1)) * cos(2.0 * 3.14159265358979323846 * u2);
  return (float)(mean + std * z);
}

/* =========================================================================
 * I/O
 * ====================================================================== */
uint8_t *read_binary_data(const char *filename, size_t num_freq,
                          size_t *out_size) {
  FILE *f = fopen(filename, "rb");
  if (!f)
    FATALF("Cannot open file: %s", filename);
  fseek(f, 0, SEEK_END);
  long file_size = ftell(f);
  fseek(f, 0, SEEK_SET);
  if (file_size < 0)
    FATAL("ftell failed");
  if ((size_t)file_size % num_freq != 0)
    FATAL("File size not multiple of num_freq");
  uint8_t *data = (uint8_t *)malloc((size_t)file_size);
  if (!data)
    FATAL("malloc failed in read_binary_data");
  if (fread(data, 1, (size_t)file_size, f) != (size_t)file_size)
    FATAL("Error reading file");
  fclose(f);
  *out_size = (size_t)file_size;
  return data;
}

void write_binary_data(const char *filename, const uint8_t *data,
                       size_t data_size, size_t nsamples, size_t nchans) {
  FILE *f = fopen(filename, "wb");
  if (!f)
    FATALF("Cannot open file for writing: %s", filename);
  // if (data_size % nsamples != 0 || data_size % nchans != 0) FATALF("Data
  // dimensions mismatch: %s", filename);
  if (fwrite(data, 1, data_size, f) != data_size)
    FATAL("Error writing file");
  fclose(f);
}

/* =========================================================================
 * Type conversions
 * ====================================================================== */
void uint8_to_float(const uint8_t *in, float *out, size_t N, float zero_off) {
/* Pure element-wise transform — no accumulation, no ordering issue. */
#pragma omp parallel for schedule(static)
  for (size_t k = 0; k < N; ++k)
    out[k] = (float)in[k] - zero_off;
}

void float_to_uint8(float *in, uint8_t *out, size_t N, float outmean,
                    float outstd) {
  int max_threads = omp_get_max_threads();
  double *t_mean = (double *)calloc((size_t)max_threads, sizeof(double));
  double *t_sq = (double *)calloc((size_t)max_threads, sizeof(double));
  if (!t_mean || !t_sq)
    FATAL("malloc failed float_to_uint8");

  int actual_threads = 1;
#pragma omp parallel shared(actual_threads)
  {
    int tid = omp_get_thread_num();
    double lmean = 0.0, lsq = 0.0;

#pragma omp for schedule(static) nowait
    for (size_t k = 0; k < N; ++k) {
      double v = in[k];
      lmean += v;
      lsq += v * v;
    }
    t_mean[tid] = lmean;
    t_sq[tid] = lsq;

#pragma omp single
    actual_threads = omp_get_num_threads();
  }

  KahanSum ks_mean = {0, 0}, ks_sq = {0, 0};
  for (int t = 0; t < actual_threads; ++t) {
    kahan_add(&ks_mean, t_mean[t]);
    kahan_add(&ks_sq, t_sq[t]);
  }
  free(t_mean);
  free(t_sq);

  double tmpmean = ks_mean.sum / (double)N;
  double tmpstd = ks_sq.sum / (double)N - tmpmean * tmpmean;
  tmpstd = sqrt(tmpstd);

  float scl = outstd / (float)tmpstd;
  float offs = outmean - scl * (float)tmpmean;

#pragma omp parallel for schedule(static)
  for (size_t k = 0; k < N; ++k) {
    float tmp = scl * in[k] + offs;
    tmp = roundf(tmp);
    tmp = clampf(tmp, 0.0f, 255.0f);
    out[k] = (uint8_t)tmp;
  }
}

/* =========================================================================
 * Comparators
 * ====================================================================== */
static int cmp_float_asc(const void *a, const void *b) {
  float fa = *(const float *)a, fb = *(const float *)b;
  return (fa > fb) - (fa < fb);
}
static int cmp_float_desc(const void *a, const void *b) {
  return -cmp_float_asc(a, b);
}
static int cmp_double_asc(const void *a, const void *b) {
  double fa = *(const double *)a, fb = *(const double *)b;
  return (fa > fb) - (fa < fb);
}
static int cmp_double_desc(const void *a, const void *b) {
  return -cmp_double_asc(a, b);
}

/* =========================================================================
 * skf_filter — Statistical Kurtosis Filter
 * ====================================================================== */
void skf_filter(float *data, float thresig, size_t nsamples, size_t nchans) {
  if (nsamples == 0 || nchans == 0)
    FATAL("Empty data in skf_filter");

  /* ---- Step 1: raw moments per channel (each channel independent) ---- */
  double *chmean1 = (double *)calloc(nchans, sizeof(double));
  double *chmean2 = (double *)calloc(nchans, sizeof(double));
  double *chmean3 = (double *)calloc(nchans, sizeof(double));
  double *chmean4 = (double *)calloc(nchans, sizeof(double));
  double *chcorr = (double *)calloc(nchans, sizeof(double));
  if (!chmean1 || !chmean2 || !chmean3 || !chmean4 || !chcorr)
    FATAL("malloc failed skf step1");

#pragma omp parallel for schedule(static)
  for (size_t j = 0; j < nchans; ++j) {
    double m1 = 0, m2 = 0, m3 = 0, m4 = 0, corr = 0;
    double last = data[j];
    for (size_t i = 0; i < nsamples; ++i) {
      double v = data[i * nchans + j];
      double v2 = v * v;
      m1 += v;
      m2 += v2;
      m3 += v2 * v;
      m4 += v2 * v2;
      corr += v * last;
      last = v;
    }
    chmean1[j] = m1;
    chmean2[j] = m2;
    chmean3[j] = m3;
    chmean4[j] = m4;
    chcorr[j] = corr;
  }

  /* ---- Step 2: derived statistics (each channel independent) ---- */
  float *chmean = (float *)calloc(nchans, sizeof(float));
  float *chstd = (float *)calloc(nchans, sizeof(float));
  float *chskewness = (float *)calloc(nchans, sizeof(float));
  float *chkurtosis = (float *)calloc(nchans, sizeof(float));
  if (!chmean || !chstd || !chskewness || !chkurtosis)
    FATAL("malloc failed skf step2");

#pragma omp parallel for schedule(static)
  for (size_t j = 0; j < nchans; ++j) {
    chmean1[j] /= nsamples;
    chmean2[j] /= nsamples;
    chmean3[j] /= nsamples;
    chmean4[j] /= nsamples;
    chcorr[j] /= (nsamples - 1);

    double mu = chmean1[j];
    double mu2 = mu * mu;
    double var = chmean2[j] - mu2;
    chmean[j] = (float)mu;

    if (var > 0.0) {
      chskewness[j] =
          (float)(chmean3[j] - 3.0 * chmean2[j] * mu + 2.0 * mu2 * mu);
      chkurtosis[j] = (float)(chmean4[j] - 4.0 * chmean3[j] * mu +
                              6.0 * chmean2[j] * mu2 - 3.0 * mu2 * mu2);
      chkurtosis[j] /= (float)(var * var);
      chkurtosis[j] -= 3.0f;
      chskewness[j] /= (float)(var * sqrt(var));
      chcorr[j] -= mu2;
      chcorr[j] /= var;
      chstd[j] = sqrtf((float)var);
    } else {
      chstd[j] = 1.0f;
      chkurtosis[j] = FLT_MAX;
      chskewness[j] = FLT_MAX;
      chcorr[j] = FLT_MAX;
    }
  }

  size_t q_idx = nchans / 4;
  float kurtosis_q1, kurtosis_q3, kurtosis_R;
  float skewness_q1, skewness_q3, skewness_R;
  double corr_q1, corr_q3, corr_R;
  {
    float *tmp = (float *)malloc(nchans * sizeof(float));
    memcpy(tmp, chkurtosis, nchans * sizeof(float));
    qsort(tmp, nchans, sizeof(float), cmp_float_asc);
    kurtosis_q1 = tmp[q_idx];
    qsort(tmp, nchans, sizeof(float), cmp_float_desc);
    kurtosis_q3 = tmp[q_idx];
    kurtosis_R = kurtosis_q3 - kurtosis_q1;
    memcpy(tmp, chskewness, nchans * sizeof(float));
    qsort(tmp, nchans, sizeof(float), cmp_float_asc);
    skewness_q1 = tmp[q_idx];
    qsort(tmp, nchans, sizeof(float), cmp_float_desc);
    skewness_q3 = tmp[q_idx];
    skewness_R = skewness_q3 - skewness_q1;
    free(tmp);
  }
  {
    double *tmp = (double *)malloc(nchans * sizeof(double));
    for (size_t j = 0; j < nchans; ++j)
      tmp[j] = chcorr[j];
    qsort(tmp, nchans, sizeof(double), cmp_double_asc);
    corr_q1 = tmp[q_idx];
    qsort(tmp, nchans, sizeof(double), cmp_double_desc);
    corr_q3 = tmp[q_idx];
    corr_R = corr_q3 - corr_q1;
    free(tmp);
  }

  int8_t *weights = (int8_t *)calloc(nchans, sizeof(int8_t));
  if (!weights)
    FATAL("malloc failed skf step4");

  long kill_count = 0;
  if (thresig >= 0) {
#pragma omp parallel for reduction(+ : kill_count) schedule(static)
    for (size_t j = 0; j < nchans; ++j) {
      if (chkurtosis[j] >= kurtosis_q1 - thresig * kurtosis_R &&
          chkurtosis[j] <= kurtosis_q3 + thresig * kurtosis_R &&
          chskewness[j] >= skewness_q1 - thresig * skewness_R &&
          chskewness[j] <= skewness_q3 + thresig * skewness_R &&
          chcorr[j] >= corr_q1 - thresig * corr_R &&
          chcorr[j] <= corr_q3 + thresig * corr_R) {
        weights[j] = 1;
      } else {
        ++kill_count;
      }
    }
  } else {
#pragma omp parallel for schedule(static)
    for (size_t j = 0; j < nchans; ++j)
      weights[j] = 1;
  }

/* ---- Step 5: normalise — each row independent ---- */
#pragma omp parallel for schedule(static)
  for (size_t i = 0; i < nsamples; ++i) {
    float *row = &data[i * nchans];
    for (size_t j = 0; j < nchans; ++j)
      row[j] = (float)weights[j] * (row[j] - chmean[j]) / chstd[j];
  }

/* ---- Step 6: bad-channel fill — thread-local RNG ---- */
#pragma omp parallel
  {
    Rng rng;
    rng_init(&rng, omp_get_thread_num());
#pragma omp for schedule(static)
    for (size_t i = 0; i < nsamples; ++i)
      for (size_t j = 0; j < nchans; ++j)
        if (weights[j] == 0)
          data[i * nchans + j] = rng_normal(&rng, 0.0f, 1.0f);
  }

  free(chmean1);
  free(chmean2);
  free(chmean3);
  free(chmean4);
  free(chcorr);
  free(chmean);
  free(chstd);
  free(chskewness);
  free(chkurtosis);
  free(weights);
}

/* =========================================================================
 * patch_filter
 * ====================================================================== */
void patch_filter(float *data, size_t nsamples, size_t nchans, int filltype) {
  int *mask = (int *)calloc(nsamples, sizeof(int));
  if (!mask)
    FATAL("malloc failed patch_filter mask");

/* ---- Step 1: detect zero-variance rows ---- */
#pragma omp parallel for schedule(static)
  for (size_t i = 0; i < nsamples; ++i) {
    double sum = 0.0, sq_sum = 0.0;
    for (size_t j = 0; j < nchans; ++j) {
      double v = data[i * nchans + j];
      sum += v;
      sq_sum += v * v;
    }
    double mean = sum / nchans;
    double var = sq_sum / nchans - mean * mean;
    if (var == 0.0) {
#pragma omp atomic write
      mask[i] = 1;
      if (i != 0) {
#pragma omp atomic write
        mask[i - 1] = 1;
      }
      if (i != nsamples - 1) {
#pragma omp atomic write
        mask[i + 1] = 1;
      }
    }
  }

  long kill_count = 0;
  for (size_t i = 0; i < nsamples; ++i)
    kill_count += mask[i];

  printf("Patch Filter: Kill rate = %.6f\n",
         (double)kill_count / (double)nsamples);

  /* ---- Step 2: mean/var of non-flagged samples ----
   * NOTE: Original C++ bug preserved — `continue` always fires here,
   * so count stays 0 and chmean_patch/chvar_patch stay zero.
   * See original conversion comment for details. */

  double *chmean_patch = (double *)calloc(nchans, sizeof(double));
  double *chvar_patch = (double *)calloc(nchans, sizeof(double));
  if (!chmean_patch || !chvar_patch)
    FATAL("malloc failed patch step2");

  long count = 0;
  for (size_t i = 0; i < nsamples; ++i) {
    if (mask[i])
      printf("This time sample is flagged: %zu\n", i);
    continue; /* faithfully reproduces original bug */

    for (size_t j = 0; j < nchans; ++j) {
      chmean_patch[j] += data[i * nchans + j];
      chvar_patch[j] += data[i * nchans + j] * data[i * nchans + j];
    }
    ++count;
  }
  if (count > 0) {
    for (size_t j = 0; j < nchans; ++j) {
      chmean_patch[j] /= count;
      chvar_patch[j] =
          chvar_patch[j] / count - chmean_patch[j] * chmean_patch[j];
    }
  }

  /* ---- Step 3: fill ---- */
  if (filltype == 0) { /* mean fill */
#pragma omp parallel for schedule(static)
    for (size_t i = 0; i < nsamples; ++i) {
      if (!mask[i])
        continue;
      for (size_t j = 0; j < nchans; ++j)
        data[i * nchans + j] = (float)chmean_patch[j];
    }
  } else { /* rand fill — thread-local RNG */
    float *ch_std = (float *)malloc(nchans * sizeof(float));
    if (!ch_std)
      FATAL("malloc failed patch rand std");
    for (size_t j = 0; j < nchans; ++j)
      ch_std[j] = (float)sqrt(chvar_patch[j] > 0.0 ? chvar_patch[j] : 0.0);

#pragma omp parallel
    {
      Rng rng;
      rng_init(&rng, omp_get_thread_num());
#pragma omp for schedule(static)
      for (size_t i = 0; i < nsamples; ++i) {
        if (!mask[i])
          continue;
        for (size_t j = 0; j < nchans; ++j)
          data[i * nchans + j] =
              rng_normal(&rng, (float)chmean_patch[j], ch_std[j]);
      }
    }
    free(ch_std);
  }

  free(mask);
  free(chmean_patch);
  free(chvar_patch);
}

/* =========================================================================
 * equalization
 * ====================================================================== */
void equalization(float *data, size_t nsamples, size_t nchans, float *chmean,
                  float *chstd) {
  memset(chmean, 0, nchans * sizeof(float));
  memset(chstd, 0, nchans * sizeof(float));

  int max_threads = omp_get_max_threads();

  /* double scratch: avoids float precision loss during accumulation */
  double *acc_mean =
      (double *)calloc((size_t)max_threads * nchans, sizeof(double));
  double *acc_sq =
      (double *)calloc((size_t)max_threads * nchans, sizeof(double));
  if (!acc_mean || !acc_sq)
    FATAL("malloc failed equalization");

  int actual_threads = 1;

#pragma omp parallel shared(actual_threads)
  {
    int tid = omp_get_thread_num();
    double *lmean = acc_mean + (size_t)tid * nchans;
    double *lsq = acc_sq + (size_t)tid * nchans;

#pragma omp for schedule(static) nowait
    for (size_t i = 0; i < nsamples; ++i) {
      const float *row = &data[i * nchans];
      for (size_t j = 0; j < nchans; ++j) {
        lmean[j] += (double)row[j];
        lsq[j] += (double)row[j] * row[j];
      }
    }

#pragma omp single
    actual_threads = omp_get_num_threads();
  }

  double *dmean = (double *)calloc(nchans, sizeof(double));
  double *dsq = (double *)calloc(nchans, sizeof(double));
  double *cmean_c =
      (double *)calloc(nchans, sizeof(double)); /* Kahan compensators */
  double *csq_c = (double *)calloc(nchans, sizeof(double));
  if (!dmean || !dsq || !cmean_c || !csq_c)
    FATAL("malloc failed equalization merge");

  for (int t = 0; t < actual_threads; ++t) {
    double *lmean = acc_mean + (size_t)t * nchans;
    double *lsq = acc_sq + (size_t)t * nchans;
    for (size_t j = 0; j < nchans; ++j) {
      /* Kahan add for dmean[j] */
      double y, tmp;
      y = lmean[j] - cmean_c[j];
      tmp = dmean[j] + y;
      cmean_c[j] = (tmp - dmean[j]) - y;
      dmean[j] = tmp;
      /* Kahan add for dsq[j] */
      y = lsq[j] - csq_c[j];
      tmp = dsq[j] + y;
      csq_c[j] = (tmp - dsq[j]) - y;
      dsq[j] = tmp;
    }
  }
  free(acc_mean);
  free(acc_sq);
  free(cmean_c);
  free(csq_c);

  for (size_t j = 0; j < nchans; ++j) {
    double m = dmean[j] / (double)nsamples;
    double s2 = dsq[j] / (double)nsamples - m * m;
    chmean[j] = (float)m;
    chstd[j] = (float)sqrt(s2 > 0.0 ? s2 : 0.0);
    if (chstd[j] == 0.0f)
      chstd[j] = 1.0f;
  }
  free(dmean);
  free(dsq);

#pragma omp parallel for schedule(static)
  for (size_t i = 0; i < nsamples; ++i) {
    float *row = &data[i * nchans];
    for (size_t j = 0; j < nchans; ++j)
      row[j] = (row[j] - chmean[j]) / chstd[j];
  }
}

typedef struct {
  float *arr;
  size_t size, cap;
  int is_max;
} Heap;

static void heap_init(Heap *h, size_t cap, int is_max) {
  h->arr = (float *)malloc(cap * sizeof(float));
  h->size = 0;
  h->cap = cap;
  h->is_max = is_max;
}
static void heap_free_h(Heap *h) {
  free(h->arr);
  h->size = 0;
}
static int heap_cmp(const Heap *h, float a, float b) {
  return h->is_max ? (a > b) : (a < b);
}
static void heap_push(Heap *h, float v) {
  size_t i = h->size++;
  h->arr[i] = v;
  while (i > 0) {
    size_t p = (i - 1) / 2;
    if (heap_cmp(h, h->arr[i], h->arr[p])) {
      float t = h->arr[i];
      h->arr[i] = h->arr[p];
      h->arr[p] = t;
      i = p;
    } else
      break;
  }
}
static float heap_top(const Heap *h) { return h->arr[0]; }
static void heap_pop(Heap *h) {
  h->arr[0] = h->arr[--h->size];
  size_t i = 0;
  for (;;) {
    size_t l = 2 * i + 1, r = 2 * i + 2, best = i;
    if (l < h->size && heap_cmp(h, h->arr[l], h->arr[best]))
      best = l;
    if (r < h->size && heap_cmp(h, h->arr[r], h->arr[best]))
      best = r;
    if (best == i)
      break;
    float t = h->arr[i];
    h->arr[i] = h->arr[best];
    h->arr[best] = t;
    i = best;
  }
}
static void heap_remove(Heap *h, float v) {
  for (size_t i = 0; i < h->size; ++i) {
    if (h->arr[i] == v) {
      h->arr[i] = h->arr[--h->size];
      while (i > 0) {
        size_t p = (i - 1) / 2;
        if (heap_cmp(h, h->arr[i], h->arr[p])) {
          float t = h->arr[i];
          h->arr[i] = h->arr[p];
          h->arr[p] = t;
          i = p;
        } else
          break;
      }
      size_t j = i;
      for (;;) {
        size_t l = 2 * j + 1, r = 2 * j + 2, best = j;
        if (l < h->size && heap_cmp(h, h->arr[l], h->arr[best]))
          best = l;
        if (r < h->size && heap_cmp(h, h->arr[r], h->arr[best]))
          best = r;
        if (best == j)
          break;
        float t = h->arr[j];
        h->arr[j] = h->arr[best];
        h->arr[best] = t;
        j = best;
      }
      return;
    }
  }
}
static void heap_rebalance(Heap *low, Heap *high) {
  if (low->size > high->size + 1) {
    heap_push(high, heap_top(low));
    heap_pop(low);
  } else if (high->size > low->size + 1) {
    heap_push(low, heap_top(high));
    heap_pop(high);
  }
}
static float heap_median(const Heap *low, const Heap *high) {
  if (low->size > high->size)
    return heap_top(low);
  if (high->size > low->size)
    return heap_top(high);
  return (heap_top(low) + heap_top(high)) / 2.0f;
}

static void sliding_median_float(const float *data, float *out, long size,
                                 int w) {
  if (w > (int)size)
    w = (int)size;
  Heap low, high;
  heap_init(&low, (size_t)w + 2, 1);
  heap_init(&high, (size_t)w + 2, 0);
  int a = -w / 2 - 1, b = (w - 1) / 2;
  heap_push(&low, data[0]);
  float median = data[0];
  for (int i = 1; i < b && i < (int)size; ++i) {
    if (data[i] >= median)
      heap_push(&high, data[i]);
    else
      heap_push(&low, data[i]);
    heap_rebalance(&low, &high);
    median = heap_median(&low, &high);
  }
  for (int i = 0; i < w / 2 + 1 && b < (int)size; ++i) {
    if (data[b] >= median)
      heap_push(&high, data[b]);
    else
      heap_push(&low, data[b]);
    heap_rebalance(&low, &high);
    median = heap_median(&low, &high);
    out[i] = median;
    ++a;
    ++b;
  }
  for (int i = w / 2 + 1; i < (int)size - (w - 1) / 2; ++i) {
    if (data[b] >= median)
      heap_push(&high, data[b]);
    else
      heap_push(&low, data[b]);
    int fl = 0;
    for (size_t k = 0; k < low.size; ++k)
      if (low.arr[k] == data[a]) {
        fl = 1;
        break;
      }
    if (fl)
      heap_remove(&low, data[a]);
    else
      heap_remove(&high, data[a]);
    heap_rebalance(&low, &high);
    median = heap_median(&low, &high);
    out[i] = median;
    ++a;
    ++b;
  }
  for (int i = (int)size - (w - 1) / 2; i < (int)size; ++i) {
    int fl = 0;
    for (size_t k = 0; k < low.size; ++k)
      if (low.arr[k] == data[a]) {
        fl = 1;
        break;
      }
    if (fl)
      heap_remove(&low, data[a]);
    else
      heap_remove(&high, data[a]);
    heap_rebalance(&low, &high);
    median = heap_median(&low, &high);
    out[i] = median;
    ++a;
  }
  heap_free_h(&low);
  heap_free_h(&high);
}

static void sliding_median_double(const double *data, double *out, long size,
                                  int w) {
  float *fi = (float *)calloc((size_t)size, sizeof(float));
  float *fo = (float *)malloc((size_t)size * sizeof(float));
  if (!fi || !fo)
    FATAL("malloc failed sliding_median_double");
  for (long i = 0; i < size; ++i)
    fi[i] = (float)data[i];
  sliding_median_float(fi, fo, size, w);
  for (long i = 0; i < size; ++i)
    out[i] = fo[i];
  free(fi);
  free(fo);
}

/* =========================================================================
 * baseline_filter
 * ====================================================================== */
void baseline_filter(float *data, size_t nsamples, size_t nchans, float width,
                     float tsamp) {
  double *s = (double *)calloc(nsamples, sizeof(double));
  if (!s)
    FATAL("malloc failed baseline s");

  int window_size = (int)(width / tsamp);
  int max_threads = omp_get_max_threads();

/* ---- Step 1: channel mean per time sample (each row independent) ---- */
#pragma omp parallel for schedule(static)
  for (size_t i = 0; i < nsamples; ++i) {
    double sum = 0.0;
    const float *row = &data[i * nchans];
    for (size_t j = 0; j < nchans; ++j)
      sum += row[j];
    s[i] = sum / nchans;
  }

  /* ---- Step 2: sliding median (inherently serial) ---- */
  double *s_med = (double *)malloc(nsamples * sizeof(double));
  if (!s_med)
    FATAL("malloc failed baseline s_med");
  sliding_median_double(s, s_med, (long)nsamples, window_size);
  memcpy(s, s_med, nsamples * sizeof(double));
  free(s_med);

  /* ---- Step 3: regression accumulation ---- */
  double *t_xe = (double *)calloc((size_t)max_threads * nchans, sizeof(double));
  double *t_xs = (double *)calloc((size_t)max_threads * nchans, sizeof(double));
  double *t_se = (double *)calloc((size_t)max_threads, sizeof(double));
  double *t_ss = (double *)calloc((size_t)max_threads, sizeof(double));
  if (!t_xe || !t_xs || !t_se || !t_ss)
    FATAL("malloc failed baseline regression");

  int actual_threads = 1;

#pragma omp parallel shared(actual_threads)
  {
    int tid = omp_get_thread_num();
    double *lxe = t_xe + (size_t)tid * nchans;
    double *lxs = t_xs + (size_t)tid * nchans;
    double lse = 0.0, lss = 0.0;

#pragma omp for schedule(static) nowait
    for (size_t i = 0; i < nsamples; ++i) {
      const float *row = &data[i * nchans];
      double si = s[i];
      for (size_t j = 0; j < nchans; ++j) {
        lxe[j] += (double)row[j];
        lxs[j] += (double)row[j] * si;
      }
      lse += si;
      lss += si * si;
    }
    t_se[tid] = lse;
    t_ss[tid] = lss;

#pragma omp single
    actual_threads = omp_get_num_threads();
  }

  double *xe = (double *)calloc(nchans, sizeof(double));
  double *xs = (double *)calloc(nchans, sizeof(double));
  double *cxe = (double *)calloc(nchans, sizeof(double)); /* compensators */
  double *cxs = (double *)calloc(nchans, sizeof(double));
  if (!xe || !xs || !cxe || !cxs)
    FATAL("malloc failed baseline merge");

  double se = 0.0, ss = 0.0, cse = 0.0, css = 0.0;
  for (int t = 0; t < actual_threads; ++t) {
    double *lxe = t_xe + (size_t)t * nchans;
    double *lxs = t_xs + (size_t)t * nchans;
    /* Kahan merge for xe, xs per channel */
    for (size_t j = 0; j < nchans; ++j) {
      double y, tmp;
      y = lxe[j] - cxe[j];
      tmp = xe[j] + y;
      cxe[j] = (tmp - xe[j]) - y;
      xe[j] = tmp;
      y = lxs[j] - cxs[j];
      tmp = xs[j] + y;
      cxs[j] = (tmp - xs[j]) - y;
      xs[j] = tmp;
    }
    /* Kahan merge for se, ss */
    {
      double y = t_se[t] - cse;
      double tmp = se + y;
      cse = (tmp - se) - y;
      se = tmp;
    }
    {
      double y = t_ss[t] - css;
      double tmp = ss + y;
      css = (tmp - ss) - y;
      ss = tmp;
    }
  }
  free(t_xe);
  free(t_xs);
  free(t_se);
  free(t_ss);
  free(cxe);
  free(cxs);

  /* ---- Step 4: coefficients (serial, cheap) ---- */
  double denom = se * se - ss * (double)nsamples;
  double *alpha = (double *)calloc(nchans, sizeof(double));
  double *beta = (double *)calloc(nchans, sizeof(double));
  if (!alpha || !beta)
    FATAL("malloc failed baseline coeffs");

  if (denom != 0.0) {
    for (size_t j = 0; j < nchans; ++j) {
      alpha[j] = (xe[j] * se - xs[j] * (double)nsamples) / denom;
      beta[j] = (xs[j] * se - xe[j] * ss) / denom;
    }
  }

/* ---- Step 5: subtract baseline (each row independent) ---- */
#pragma omp parallel for schedule(static)
  for (size_t i = 0; i < nsamples; ++i) {
    float *row = &data[i * nchans];
    double si = s[i];
    for (size_t j = 0; j < nchans; ++j)
      row[j] -= (float)(alpha[j] * si + beta[j]);
  }

  free(s);
  free(xe);
  free(xs);
  free(alpha);
  free(beta);
}

/* =========================================================================
 * process_block
 * ====================================================================== */
static void process_block(uint8_t *block_ptr, size_t block_N, size_t nsamples,
                          size_t nchans, int filltype, float thresig,
                          float outmean, float outstd, float width,
                          float tsamp) {
  float *fd = (float *)malloc(block_N * sizeof(float));
  float *cmean = (float *)malloc(nchans * sizeof(float));
  float *cstd = (float *)malloc(nchans * sizeof(float));
  if (!fd || !cmean || !cstd)
    FATAL("malloc failed in process_block");

  uint8_to_float(block_ptr, fd, block_N, -64.0f);
  // patch_filter(fd, nsamples, nchans, filltype);
  skf_filter(fd, thresig, nsamples, nchans);
  equalization(fd, nsamples, nchans, cmean, cstd);
  baseline_filter(fd, nsamples, nchans, width, tsamp);
  float_to_uint8(fd, block_ptr, block_N, outmean, outstd);

  free(fd);
  free(cmean);
  free(cstd);
}

int initialize_TEL_SHM(void) {
  int id_TEL_hdr = 0, id_TEL_buf = 0;
  int CurrSamples;
  long int CurrTotalWords;
  int output_size = 2; // 2 for 8bit data, hardcoded.

  // bad chans ranges = 0-200     : 550 - 559.765625
  //                  = 3000-4095 : 696.484375 - 750;

  /*
          int nbadchRange_B3;
          int nbadchRange_B4;
          int nbadchRange_B5;

          nbadchRange_B4=2;

          B3_Chan_range[10];
          B4_Chan_range[10];
          B5_Chan_range[10];
  */

  id_TEL_hdr = shmget(BeamHdrKey, sizeof(BeamHeaderType), SHM_RDONLY);
  if (id_TEL_hdr < 0) {
    perror("TEL shmget BEAM HEADER ID");
    return -1;
  }

  dataHdr_TEL = (BeamHeaderType *)shmat(id_TEL_hdr, 0, SHM_RDONLY);
  if (dataHdr_TEL == (void *)-1) {
    perror("TEL shmat BEAM HEADER");
    return -1;
  }

  nBeams = dataHdr_TEL->BeamGenHdr.BeamSteeringParams.nSteeringBeamsPerNode;
  // nBeams = 1; // nBeams hardcoded to be 1 for IA beam hosts.
  curr_chans = dataHdr_TEL->corr.corrpar.channels;
  Freq = dataHdr_TEL->ScanTab[0].source.freq[0];

  if (myrank == 0) {
    fprintf(stderr, "##########################################################"
                    "################\n");
    fprintf(stderr, "RF = %lf\n", Freq);
    fprintf(stderr, "##########################################################"
                    "################\n");
  }
  /*
          ChWidth  = dataHdr_TEL->ScanTab[0].source.ch_width * -1 *
     dataHdr_TEL->ScanTab[0].source.net_sign[0];

          for (int ch=0;ch<curr_chans;ch++) {
                  chan_freq = Freq + ChWidth*ch;
          }
  */

  if (myrank == 0)
    fprintf(stderr, "############# CHANS = %d \n", curr_chans);

  CurrSamples = (32 * 25); // 800 (no. of FFT in 1Sec. for 4096 channel)
  // CurrSamples = (32 * 100);				 // 3200 (no. of FFT in
  // 1Sec. for 1024 channel)
  CurrTotalWords = CurrSamples * curr_chans; // independent of current channels.
  CurrRecSize =
      CurrTotalWords * WordSize / (output_size); // WordSie=sizeof(short)
  // CurrRecSize     = CurrTotalWords*WordSize;
  // CurrRecSize     *= nBeams;
  dataSize =
      TEL_to_FRB_Block_factor * CurrRecSize; // FRB SHM size of single beam.
  // dataSize = CurrRecSize; // FRB SHM size of single beam.
  if (myrank == 0)
    fprintf(stderr,
            "############# nBeams = %d CurrRecSize = %ld dataSize = %ld \n",
            nBeams, CurrRecSize, dataSize);

  int local_ExtraBuf = 64;
  // long CurrShmDataSize = (((MaxRecs + 1) * TEL_SHM_BlockSize * nBeams +
  // local_ExtraBuf) / PageSize + 1) * PageSize;
  long CurrShmDataSize =
      (((MaxRecs + 1) * CurrRecSize * nBeams + local_ExtraBuf) / PageSize + 1) *
      PageSize;
  int local_ShmDataOff =
      ((sizeof(GlobalInfoType) + local_ExtraBuf) / PageSize + 1) * PageSize;
  long CurrShmSize = local_ShmDataOff + CurrShmDataSize;

  id_TEL_buf = shmget(ShmKey, CurrShmSize, SHM_RDONLY);
  if (id_TEL_buf < 0) {
    perror("TEL shmget BEAM BUFFER ID");
    return -1;
  }
  dataBuffer_TEL = (GlobalInfoType *)shmat(id_TEL_buf, 0, SHM_RDONLY);
  if ((void *)dataBuffer_TEL == (void *)-1) {
    perror("TEL shmat BEAM BUFFER ID");
    return -3;
  }
  Shmp = (unsigned char *)dataBuffer_TEL;

  Rec = dataBuffer_TEL->rec;
  RecNum = ((unsigned int)(dataBuffer_TEL->rec_ind - 1)) % MaxRecs;
  // RecNum = ((unsigned int) (dataBuffer_TEL->rec_ind - MaxRecs + 1) %
  // MaxRecs);
  if (myrank == 0)
    fprintf(stderr, "dataBuffer_TEL->rec_ind: %d RecNum = %d \n",
            dataBuffer_TEL->rec_ind, RecNum);
  // while(RecNum != 2 ) RecNum = (((unsigned int) (dataBuffer_TEL->rec_ind -
  // 1)) % MaxRecs); fprintf (stderr, "RecNum: %u\n", RecNum);
  DataSeq = Rec[RecNum].rec_seq; /* expected first rec */

  return 0;
}

int initialize_FRB_SHM() {
  int id_FRB_hdr, id_FRB_buf;
  int curr_chans_FRB, nBeams_FRB;

  id_FRB_hdr = shmget(DasHeaderKey, sizeof(BeamHeaderType), IPC_CREAT | 0666);
  if (id_FRB_hdr < 0) {
    perror("FRB shmget HEADER ID");
    return -1;
  }

  dataHdr_FRB = (BeamHeaderType *)shmat(id_FRB_hdr, 0, 0);
  if ((void *)dataHdr_FRB == (void *)-1) {
    perror("FRB shmat HEADER ID");
    return -3;
  }

  // nBeams_FRB =
  // dataHdr_FRB->BeamGenHdr.BeamSteeringParams.nSteeringBeamsPerNode;
  // curr_chans_FRB = dataHdr_FRB->corr.corrpar.channels;

  // CurrSamples = (32*25);  // 800 (no. of FFT in 1Sec.)
  // CurrTotalWords = CurrSamples * curr_chans_FRB  ; // independent of current
  // channels. CurrRecSize     = CurrTotalWords*WordSize/(output_size); //
  // WordSie=sizeof(short) CurrRecSize     = CurrTotalWords*WordSize;
  // CurrRecSize     *= nBeams_FRB;
  // dataSize = TEL_to_FRB_Block_factor * CurrRecSize;
  // FRB_SHM_Size = dataSize*MaxDataBlocks;

  id_FRB_buf = shmget(DasBufferKey, sizeof(DataBuffer), IPC_CREAT | 0666);
  if (id_FRB_buf < 0) {
    fprintf(stderr, "DATABUFFER SIZE = %ld\n", (long)sizeof(DataBuffer));
    perror("FRB shmget BUFFER Problem :");
    return -1;
  }
  dataBuffer_FRB = (DataBuffer *)shmat(id_FRB_buf, 0, 0);
  if ((void *)dataBuffer_FRB == (void *)-1) {
    perror("FRB shmat BUFFER ID");
    return -3;
  }

  // Copying the TEL_SHM header as it is for right now. This needs to be updated
  // multiple times during the observation.
  memcpy(dataHdr_FRB, dataHdr_TEL, sizeof(BeamHeaderType));
  // UP: Need to make both total and per node of beams to 1, since both AA and
  // shazam rely on these variables to calculate SHM size, and connect to SHM.
  // dataHdr_FRB->BeamGenHdr.BeamSteeringParams.nSteeringBeams = 1; // provision
  // for IA beam, single beam per host.
  // dataHdr_FRB->BeamGenHdr.BeamSteeringParams.nSteeringBeamsPerNode = 1; //
  // provision for IA beam, single beam per host.

  bandsign = dataHdr_TEL->ScanTab[0].source.net_sign[0];
  if (bandsign == -1) { // i.e., net_sign=-1, means USB (we need band flip)
    dataHdr_FRB->ScanTab[0].source.net_sign[0] *= -1;
    dataHdr_FRB->ScanTab[0].source.net_sign[1] *= -1;
    dataHdr_FRB->ScanTab[0].source.freq[0] =
        dataHdr_FRB->ScanTab[0].source.freq[0] +
        (1e6 * dataHdr_FRB->corr.daspar.gsb_acq_bw) /
            dataHdr_FRB->corr.daspar.gsb_final_bw;
    dataHdr_FRB->ScanTab[0].source.freq[1] =
        dataHdr_FRB->ScanTab[0].source.freq[1] +
        (1e6 * dataHdr_FRB->corr.daspar.gsb_acq_bw) /
            dataHdr_FRB->corr.daspar.gsb_final_bw;
    // if (myrank == 0) fprintf(stderr, "#### BAND SIGN = %d freq1 = %lf, freq2
    // = %lf\n", bandsign,
    // dataHdr_FRB->ScanTab[0].source.freq[0],dataHdr_FRB->ScanTab[0].source.freq[1]);
    //		dataHdr_FRB->ScanTab[0].source.ch_width = -1 *
    //dataHdr_FRB->ScanTab[0].source.ch_width; // already header is copied,
    //which has inverted sign.
  }

  bandwidth = dataHdr_FRB->corr.daspar.gsb_acq_bw /
              dataHdr_FRB->corr.daspar.gsb_final_bw;

  dataBuffer_FRB->timestamp_gps[0] = Rec[RecNum].timestamp_gps;
  dataBuffer_FRB->blk_nano[0] =
      Rec[RecNum].blk_nano; // Adding nanoseconds accuracy for the timestamps.
  dataBuffer_FRB->pc_time = Rec[RecNum].pc_time;
  dataBuffer_FRB->ref_time = dataBuffer_TEL->ref_time;
  dataBuffer_FRB->rec_time = Rec[RecNum].rec_time;

  dataBuffer_FRB->nBeams = nBeams;
  dataBuffer_FRB->curRecord = 0;
  dataBuffer_FRB->curBlock = 0;
  dataBuffer_FRB->blockSize = CurrRecSize;
  dataBuffer_FRB->active = 1;
  dataBuffer_FRB->is_buf_empty =
      1; // There is no meaningful data in the buffer. After filling in one
         // block of buffer, this flag will be set to 0.
  dataBuffer_FRB->flag = 0;

  return 0;
}

int process_acq(void) {
  if (dataBuffer_TEL->acq_flag & AcqOver) {
    fprintf(stderr, "ACQOVER\n");
    return -2;
  }

  ind = dataBuffer_TEL->blk_ind > 0 ? dataBuffer_TEL->blk_ind - 1 : MaxBLK - 1;
  if (dataBuffer_TEL->blk_seq[ind] > DataSeq + MaxRecs - 1) {
    fprintf(stderr,
            "Processing lagged behind: DataSeq=%d AcqSeq=%d, RecNum=%d.\n",
            DataSeq, dataBuffer_TEL->blk_seq[ind], RecNum);
    fprintf(stderr, "New Seq. for processing=%d\a\n",
            dataBuffer_TEL->blk_seq[ind] - MaxRecs + 1);
    DataSeq = dataBuffer_TEL->blk_seq[ind] - MaxRecs + 1;
    RecNum = (dataBuffer_TEL->rec_ind - MaxRecs + 1) % MaxRecs;
  }

  while (1) {
    if (Rec[RecNum].rec_flag & Marked) {
      return -3;
    }
    if (Rec[RecNum].rec_seq < DataSeq) {
      RecNum = (RecNum + 1) % MaxRecs;
    } else {
      break;
    }
  }

  if (Rec[RecNum].rec_seq > DataSeq) {
    DataSeq++;
    return -4;
  }

  if ((Rec[RecNum].rec_flag & GoodData) == 0) {
    fprintf(stderr, "Bad data: DataSeq=%d flag=%4d\n", DataSeq,
            Rec[RecNum].rec_flag);
    ++DataSeq;
    return -5;
  }

  blkp_TEL = Shmp + Rec[RecNum].beg_off;
  Data_Time[RecNum] = Rec[RecNum].pc_time;
  local_t = localtime(&Rec[RecNum].timestamp_gps.tv_sec);
  gettimeofday(&tv1, &tz1);
  Data_Time[RecNum] =
      ((tv1.tv_sec - tz1.tz_minuteswest * 60) % 86400 + tv1.tv_usec / 1e6);
  // Data_Time[RecNum] = ((tv1.tv_sec - 19800) % 86400 + tv1.tv_usec / 1e6);

  memcpy(dataHdr_FRB, dataHdr_TEL, sizeof(BeamHeaderType));
  // dataHdr_FRB->BeamGenHdr.BeamSteeringParams.nSteeringBeams = 1; // provision
  // for IA beam, single beam per host.
  // dataHdr_FRB->BeamGenHdr.BeamSteeringParams.nSteeringBeamsPerNode = 1; //
  // provision for IA beam, single beam per host.
  bandsign = dataHdr_TEL->ScanTab[0].source.net_sign[0];
  if (bandsign == -1) { // i.e., net_sign=-1, means USB (we need band flip)
    dataHdr_FRB->ScanTab[0].source.net_sign[0] *= -1;
    dataHdr_FRB->ScanTab[0].source.net_sign[1] *= -1;
    dataHdr_FRB->ScanTab[0].source.freq[0] =
        dataHdr_FRB->ScanTab[0].source.freq[0] +
        (1e6 * dataHdr_FRB->corr.daspar.gsb_acq_bw) /
            dataHdr_FRB->corr.daspar.gsb_final_bw;
    dataHdr_FRB->ScanTab[0].source.freq[1] =
        dataHdr_FRB->ScanTab[0].source.freq[1] +
        (1e6 * dataHdr_FRB->corr.daspar.gsb_acq_bw) /
            dataHdr_FRB->corr.daspar.gsb_final_bw;
    // if (myrank == 0) fprintf(stderr, "#### BAND SIGN = %d freq1 = %lf, freq2
    // = %lf\n", bandsign,
    // dataHdr_FRB->ScanTab[0].source.freq[0],dataHdr_FRB->ScanTab[0].source.freq[1]);
    //		dataHdr_FRB->ScanTab[0].source.ch_width = -1 *
    //dataHdr_FRB->ScanTab[0].source.ch_width; // already header is copied,
    //which has inverted sign.
  }

  if (TEL_to_FRB_RecNum == 0 &&
      dataBuffer_FRB->flag != 0) // Data for 1 block of FRB_SHM has accumulated.
                                 // Update the curRecord and curBlock.
  {
    DataSeq_FRB++;
    RecNum_FRB = (RecNum_FRB + 1) % MaxDataBlocks;

    dataBuffer_FRB->timestamp_gps[RecNum_FRB] = Rec[RecNum].timestamp_gps;
    dataBuffer_FRB->blk_nano[RecNum_FRB] =
        Rec[RecNum].blk_nano; // Adding nanoseconds accuracy for the timestamps.
    dataBuffer_FRB->pc_time = Rec[RecNum].pc_time;
    dataBuffer_FRB->ref_time = dataBuffer_TEL->ref_time;
    dataBuffer_FRB->rec_time = Rec[RecNum].rec_time;

    dataBuffer_FRB->curRecord = RecNum_FRB;
    dataBuffer_FRB->curBlock = DataSeq_FRB;
    dataBuffer_FRB->is_buf_empty = 0; // Now it is sure that there is at least
                                      // one block of buffer filled with data.

    if (myrank == 0)
      fprintf(stderr, "#--------------------------------------------\n");
    // if (myrank == 0) fprintf(stderr, "##-- Block number in FRB_SHM: %d in
    // TEL_SHM: %d\n", dataBuffer_FRB->curRecord,
    // dataBuffer_FRB->curRecord*TEL_to_FRB_Block_factor);
    if (myrank == 0)
      fprintf(stderr, "##-- Block number in FRB_SHM: %d\n",
              dataBuffer_FRB->curRecord);
    if (myrank == 0)
      fprintf(stderr, "##-- Total blocks read: %d\n", dataBuffer_FRB->curBlock);
    if (myrank == 0)
      fprintf(stderr, "#--------------------------------------------\n");
  }

  // Each block of TEL_SHM will be copied as soon as it is available in the SHM.
  // The curRecord and other parameters are updated (see the above if condition)
  // only after accumulating 32 blocks of TEL_SHM. Pretending as if FRB_SHM is
  // 12 * 32 blocks long, with each block equal in size to that of TEL_SHM.

  dataSize = TEL_to_FRB_Block_factor * CurrRecSize;
  int8_t spec_chan;
  // unsigned char tempbuffer[curr_chans];
  int8_t *spec_ptr;
  unsigned char blockbuffer[CurrRecSize];
  struct timeval start, end;
  long seconds;
  long microseconds;
  double elapsed_time;

  int fft_no, chan_no;
  size_t blk;

  gettimeofday(&start, NULL);
  // uint8_t *temp_block = malloc(nsamples * nchans * sizeof(uint8_t));
  /* FIX 2+3: Pre-allocate a contiguous buffer pool for all beams outside the
   * parallel loop.  Each thread gets its own slice beam_bufs[blk * block_N],
   * eliminating per-iteration malloc/free and the associated allocator-lock
   * contention (~3.2 MB × 40 allocations per call). */
  size_t block_N = nsamples * (size_t)nchans;
  uint8_t *beam_bufs =
      (uint8_t *)malloc((size_t)nBeams * block_N * sizeof(uint8_t));
  if (!beam_bufs)
    FATAL("malloc failed for beam_bufs");

  // ################################################################

  /* FIX 1: blkp_FRB was a global pointer written by every thread → data race.
   *         It is now a local pointer inside the loop body.
   * FIX 2: local_block is a slice of the pre-allocated pool — no malloc/free.
   * FIX 3 (nested OMP): process_block() contains its own omp parallel regions.
   *         Running them inside an already-active 40-thread team causes up to
   *         40×40 = 1600 threads competing for cores.  We disable nested
   *         parallelism so the inner parallel regions become no-ops while the
   *         outer beam loop drives parallelism. */
  omp_set_nested(0); /* inner omp regions run single-threaded */
  omp_set_max_active_levels(1);

#pragma omp parallel for num_threads(16) schedule(static) shared(blkp_TEL)     \
    private(spec_ptr, spec_chan, chan_no, fft_no)
  for (beam_ID = 0; beam_ID < (size_t)nBeams; beam_ID++) {
    // printf("Threads running = %d\n", omp_get_num_threads());
    /* Each beam has its own private output pointer — no race condition. */
    unsigned char *local_blkp_FRB = dataBuffer_FRB->data +
                                    ((long)(dataSize * nBeams * RecNum_FRB)) +
                                    ((long)dataSize * (long)beam_ID) +
                                    ((long)CurrRecSize * TEL_to_FRB_RecNum);

    /* Each thread works on its own non-overlapping slice of beam_bufs. */
    uint8_t *local_block = beam_bufs + beam_ID * block_N;

    // fprintf(stderr, "BLK = %zu Thread %d of %d\n", beam_ID,
    // omp_get_thread_num(), omp_get_num_threads());
    memcpy(local_block, blkp_TEL + beam_ID * block_N,
           block_N * sizeof(uint8_t));

    double blk_start = omp_get_wtime() * 1e6;

    process_block(local_block, block_N, nsamples, (size_t)nchans, filltype,
                  thresig, outmean, outstd, width, tsamp);

    double blk_time = omp_get_wtime() * 1e6 - blk_start;
    float rtf = (block_size * tsamp * 1e6f) / (float)blk_time;
    (void)rtf;            /* suppress unused-variable warning; log if needed */
                          // #########################################
    if (bandsign == -1) { // net_sign=1, means LSB (we assume no band flip)
                          // fprintf(stderr, "BND INV applied...\n");
      for (fft_no = 0; fft_no < (32 * 25); fft_no++) {
        spec_ptr = local_block + fft_no * curr_chans;
        for (chan_no = 0; chan_no < (int)(curr_chans / 2); chan_no++) {
          spec_chan = *(spec_ptr + chan_no);
          *(spec_ptr + chan_no) = *(spec_ptr + curr_chans - chan_no - 1);
          *(spec_ptr + (curr_chans - chan_no - 1)) = spec_chan;
          // spec_chan = spec_ptr[chan_no];
          // spec_ptr[chan_no] = spec_ptr[curr_chans - chan_no - 1];
          // spec_ptr[curr_chans - chan_no - 1] = spec_chan;
        }
      }
    }
    // #########################################
    memcpy(local_blkp_FRB, local_block, (size_t)CurrRecSize);
  }

  free(beam_bufs);

  /*
          for (beam_ID = 0; beam_ID < nBeams; beam_ID++) {
             blkp_FRB = dataBuffer_FRB->data + ((long)(dataSize * nBeams *
     RecNum_FRB)) + ((long) dataSize * beam_ID) + ((long)CurrRecSize *
     TEL_to_FRB_RecNum); for (fft_no = 0; fft_no < (32 * 25); fft_no++) {
                          memcpy(tempbuffer, block_ptr + fft_no * curr_chans,
     curr_chans);
                          //memcpy(tempbuffer, block_ptr + (long)CurrRecSize *
     beam_ID + fft_no * curr_chans, curr_chans); if (bandsign == -1) { //
     net_sign=1, means LSB (we assume no band flip)
                                  //if (myrank == 0) fprintf(stderr, "Fliping
     FREQ = %lf -- ", dataHdr_FRB->ScanTab[0].source.freq[0]); for (chan_no = 0;
     chan_no < (int)(curr_chans / 2); chan_no++) { temp = tempbuffer[chan_no];
                                          tempbuffer[chan_no] =
     tempbuffer[curr_chans - chan_no - 1]; tempbuffer[curr_chans - chan_no - 1]
     = temp;
                                  }
                          }
                          for (chan_no = 0; chan_no < curr_chans; chan_no++) {
     blockbuffer[chan_no + fft_no * curr_chans] = (unsigned char)((short
     int)tempbuffer[chan_no]); }
          }
          memcpy(blkp_FRB, blockbuffer, CurrRecSize);
          //memcpy(blkp_FRB, block_ptr, CurrRecSize);
          //memcpy(blkp_FRB, blkp_TEL + (long)CurrRecSize * beam_ID,
     CurrRecSize);
       }
  */
  gettimeofday(&end, NULL);
  seconds = end.tv_sec - start.tv_sec;
  microseconds = end.tv_usec - start.tv_usec;
  elapsed_time = seconds + microseconds / 1000000.0;

  if (myrank == 0)
    fprintf(stderr, "## -- FILT+BAND FLIP : %lf seconds. : ", elapsed_time);
  usec = (Rec[RecNum].timestamp_gps.tv_usec);
  blk_nano = (Rec[RecNum].blk_nano);
  strftime(time_string, sizeof(time_string), "%Y %m %d %H %M %S", local_t);
  if (myrank == 0)
    fprintf(stderr,
            "##-- Current data timestamp = %02d:%02d:%02d.%06d%03d \t "
            "Rec[RecNum].rec_seq: %u \t RecNum: %u\n",
            local_t->tm_hour, local_t->tm_min, local_t->tm_sec, usec,
            (int)(1 * blk_nano), Rec[RecNum].rec_seq, RecNum);

  DataSeq++;
  RecNum = (RecNum + 1) % MaxRecs;
  TEL_to_FRB_RecNum =
      (TEL_to_FRB_RecNum + 1) %
      TEL_to_FRB_Block_factor; // TEL_to_FRB_Block_factor=32 (frb_shm.h)
  dataBuffer_FRB->flag = 1;
  return 0;
}

int main(int argc, char *argv[]) {
  struct timespec req;
  long obs_duration = round(3600 / 1.048576);

  if (argc != 2) {
    fprintf(stderr, "\nUSAGE: %s [observation duration (seconds) ]\n\n",
            argv[0]);
    fprintf(stderr, "Reading 1 hour of data by default.\n\n");
    return (-1);
  } else
    obs_duration = round(atoi(argv[1]) / 1.048576);

  req.tv_sec = 0;
  req.tv_nsec = 100000000L;

  block_size = 800;
  thresig = 4.0;
  filltype = 1;
  nchans = 4096;
  outmean = 64.0f;
  outstd = 3.0f;
  width = 0.5f;
  tsamp = 1.31072e-3f;
  raw_size = nBeams * 800 * 4096;
  nsamples = (size_t)block_size;
  block_len = nsamples * (size_t)nchans;
  n_full = raw_size / block_len; // this will always be 1

  // printf("Total blocks: %zu  |  Remainder bytes: %zu\n", n_full, remainder);
  // printf("Total blocks: %zu  \n", n_full);

  /* FIX 6: Use MPI_Init_thread so MPI is aware of our OpenMP usage.
   * MPI_THREAD_FUNNELED means only the master thread makes MPI calls,
   * which matches how this code uses MPI (Barrier in the main loop only).
   * omp_set_num_threads is moved to AFTER MPI_Init so that MPI's internal
   * thread-affinity setup does not interfere with our thread count. */
  int mpi_thread_provided;
  MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &mpi_thread_provided);
  MPI_Comm_size(MPI_COMM_WORLD, &commSize);
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

  /* Set OMP threads AFTER MPI init — safe for all MPI implementations.
   * Nested parallelism is disabled globally; inner omp regions in
   * process_block() will run single-threaded while the outer beam loop
   * drives the 40-way parallelism. */
  omp_set_num_threads(16);
  omp_set_nested(0);
  omp_set_max_active_levels(1);
  if (myrank == 0)
    fprintf(stderr, "OpenMP threads: %d  (nested parallelism disabled)\n",
            omp_get_max_threads());

  if (initialize_TEL_SHM() < 0) {
    fprintf(stderr, "TEL_SHM initialization failure\n");
    return -3;
  }
  if (initialize_FRB_SHM() < 0) {
    fprintf(stderr, "FRB_SHM initialization failure\n");
    return -3;
  }

  if (myrank == 0)
    fprintf(stderr, "----------------------------------------------------------"
                    "---------------------------\n");
  if (myrank == 0)
    fprintf(stderr,
            "   Num. Beams = %d TEL_SHM BLOCK SIZE = %ld, FRB_SHM BLOCK SIZE = "
            "%ld\n",
            (int)nBeams, (long)CurrRecSize, (long)(nBeams * (long)dataSize));
  if (myrank == 0)
    fprintf(stderr, "----------------------------------------------------------"
                    "---------------------------\n");

  int done;
  long loop_cntr = 0;

  block_ptr = (uint8_t *)malloc(nsamples * nchans * sizeof(uint8_t));

  while (1) {
    done = process_acq();
    // if (done < 0) { usleep(100000); }
    if (done < 0) {
      nanosleep(&req, NULL);
    } else if (loop_cntr++ >= obs_duration)
      return 0;
    MPI_Barrier(MPI_COMM_WORLD);
  }
  MPI_Finalize();
  return 0;
}