DP3 bandwidth / time smearing test on DUG

dug_smearing.md

This didn't work but I'm documenting it here in case I want to try it again.

take calibration solutions, fit linear model, update cable lengths in metafits before a second round of preprocessing.

singularity shell --cleanenv --bind /data/curtin_mwaeor/src/mwa-demo/demo:/demo --bind $PWD --workdir $PWD /data/curtin_mwaeor/singularity/mwa-demo_autos_cuda12.5.1.sif
/demo/82_calfit.py \
    --metafits /data/curtin_mwaeor/dev/1274010302/raw/1274010302.metafits \
    --solns /data/curtin_mwaeor/dev/1274010302/hyp_soln_1274010302_kr_75-1667l_src1k_500it_sdc50_el36.692.fits \
    --name hyp_soln --out-dir /data/curtin_mwaeor/dev/1274010302/ \
    --plot-residual
/demo/82_update_cable_lengths.py \
    --metafits /data/curtin_mwaeor/dev/1274010302/raw/1274010302.metafits \
    --name hyp_soln \
    --out-dir /data/curtin_mwaeor/dev/1274010302/ \
    --tsv /data/curtin_mwaeor/dev/1274010302/1274010302\ hyp_soln_t000_phase_fits.tsv

# creates /data/curtin_mwaeor/dev/1274010302/raw/1274010302_fixed.metafits

preprocess

birli \
    -m /data/curtin_mwaeor/dev/1274010302/raw/1274010302_fixed.metafits \
    --avg-freq-res 80 \
    --avg-time-res 8 \
    -M /data/curtin_mwaeor/dev/1274010302/prep/1274010302_2s_40kHz.ms \
    -- /data/curtin_mwaeor/dev/1274010302/raw/1274010302_2*.fits

calibration skymodel for dp3

cd /home/dev/src/src-workloads/tasks/karabo

# 1. Ensure metafits and Global sky model exists
[[ -f data/GGSM_updated.fits ]] || make get-data MODEL=GGSM_updated.fits STEP=sim
export OBSID=1274010302
[[ -f data/$OBSID.metafits ]] || wget -q -O data/$OBSID.metafits "http://ws.mwatelescope.org/metadata/fits?obs_id=$OBSID"

# 2. Generate local sourcelist from GGSM_updated.fits
# This creates both AO and BBS format sourcelists
make TELESCOPE=MWA OBSID=1274010302 DATA=1274010302 MODEL_COUNT=2000 reduce-model-mwa-21cm-forest # EXECUTOR=singularity

# 3. copy them to cal
# mkdir -p /data/curtin_mwaeor/dev/1274010302/cal/
cp data/1274010302_2000.ao.txt data/1274010302_2000_ddecal.skymodel.txt /data/curtin_mwaeor/dev/1274010302/cal/
cp /data/curtin_mwaeor/src/src-workloads/tasks/karabo/data/*.reg /data/curtin_mwaeor/dev/1274010302/cal/

one direciton hack for sky model, to use smoothnessconstraint we need ddecal, but with a single direction

awk 'BEGIN{FS=OFS=", "} NR<=2 || NF<2 {print; next} {$2="onedirection"; print}' \
    /data/curtin_mwaeor/dev/1274010302/cal/1274010302_2000_ddecal.skymodel.txt \
    > /data/curtin_mwaeor/dev/1274010302/cal/1274010302_2000_onedirection.skymodel.txt

dp3 predict (with and without smearing corrections)

cd /data/curtin_mwaeor/dev/1274010302/
module load singularity/default
MS_IN=prep/1274010302_2s_40kHz.ms
MS_MODEL=cal/1274010302_2s_40kHz_dp3nosmear_model.ms
COEFF_H5=/data/curtin_mwaeor/mwa_full_embedded_element_pattern.h5
CAL_H5=cal/1274010302_2s_40kHz_dp3nosmear.h5
MODEL_TXT=cal/1274010302_2000_onedirection.skymodel.txt

tmpdir=$(mktemp -d /dev/shm/deleteme.XXXXXX)
MS_TMP=${tmpdir}/${MS_IN##*/}
COEFF_TMP=${tmpdir}/${COEFF_H5##*/}
cp -R "${MS_IN}" "${MS_TMP}"
cp "${COEFF_H5}" "${COEFF_TMP}"
numthreads=$(( $(nproc) / 2 ))

# Apparently you can't do it all in one pass?
#     msin=${MS_TMP} \
#     msin.datacolumn=DATA \
#     msout=${MS_IN} \
#     steps=[ddecal,applycal] \
#     ddecal.sourcedb=${MODEL_TXT} \
#     ddecal.h5parm=${CAL_H5} \
#     ddecal.keepmodel=True \
#     ddecal.applysolution=false \
#     applycal.parmdb=${CAL_H5} \
#     applycal.correction=amplitude000,phase000

/usr/bin/time -v singularity exec \
    --env OMP_NUM_THREADS=${numthreads} \
    --env OPENBLAS_NUM_THREADS=1 \
    -B/data \
    -B$tmpdir \
    -W$PWD \
    docker://d3vnull0/dp3-mwa:latest DP3 \
    numthreads=${numthreads} \
    msin=${MS_TMP} \
    msout=${MS_MODEL} \
    msout.datacolumn=DATA \
    msout.overwrite=true \
    steps=[ddecal] \
    ddecal.sourcedb=${MODEL_TXT} \
    ddecal.h5parm=${CAL_H5} \
    ddecal.usebeammodel=true \
    ddecal.uvlambdamin=75 \
    ddecal.uvlambdamax=1667 \
    ddecal.solint=20 \
    ddecal.smoothnessconstraint=4e6 \
    ddecal.beamproximitylimit=300 \
    ddecal.maxiter=500 \
    ddecal.tolerance=1e-20 \
    ddecal.coefficients_path=${COEFF_TMP} \
    ddecal.keepmodel=True
    # ddecal.correctfreqsmearing=True \
    # ddecal.correcttimesmearing=True \

# apply solutions
/usr/bin/time -v singularity exec \
    --env OMP_NUM_THREADS=${numthreads} \
    --env OPENBLAS_NUM_THREADS=1 \
    -B/data -B$tmpdir -W$PWD \
    docker://d3vnull0/dp3-mwa:latest DP3 \
    numthreads=${numthreads} \
    msin=${MS_TMP} \
    msout=${MS_IN} \
    msout.datacolumn=CORRECTED_DATA \
    msout.overwrite=true \
    steps=[applycal] \
    applycal.parmdb=${CAL_H5} \
    applycal.correction=amplitude000,phase000


# plot solutions

# Create plots with losoto (amplitude and phase vs freq, per antenna)
mkdir -p plots && cat > losoto-bp.parset << 'EOF'
Ncpu = 0

[plotA_freq]
operation = PLOT
soltab = sol000/amplitude000
time.minmaxstep = [0,1e20,1000]
axesInPlot = freq
axisInTable = ant
axisInCol = pol
plotFlag = True
prefix = plots/amp_
minmax = [0,2]

[plotP_freq]
operation = PLOT
soltab = sol000/phase000
time.minmaxstep = [0,1e20,1000]
axesInPlot = freq
axisInTable = ant
axisInCol = pol
plotFlag = True
prefix = plots/phase_
minmax = [-3.15,3.15]
EOF
singularity exec -B /data docker://d3vnull0/karabo-workloads:sha-967aa66-cpu /opt/view/bin/losoto ${CAL_H5} losoto-bp.parset

# Dump text summary of amplitude and phase solutions
singularity exec -B /data docker://d3vnull0/karabo-workloads:sha-967aa66-cpu python3 /scripts/dump-h5parm-to-text.py ${CAL_H5} --outdir ./text --mode summary --soltabs amplitude000 phase000

image

module load wsclean/3.5
run_wsclean() {
    chg="$1" # e.g. "14h19m08.52539063s 6d28m32.74452209s"
    ms="$2"       # ms file input, e.g. model_1274010302_sdc${sdc}_el${el}_40kHz_4s.ms or other
    column="$3"
    target="$4"          # for log file naming, e.g. "model" or "data"
    datacolumn="DATA"
    msbase=${ms##*/}
    mshead=${msbase%.ms}
    tempdir=$(mktemp -d /dev/shm/deleteme.XXXXXX)
    cp -R "${ms}" "${tempdir}/"
    chgcentre "${tempdir}/${msbase}" $chg
    nch=96
    wsclean \
        -mgain 0.8 \
        -nmiter 10 \
        -niter 100000 \
        -size 128 128 \
        -pol XXYY \
        -data-column "${datacolumn}" \
        -name "${tempdir}/${mshead}_${datacolumn}_${target}" \
        -scale 10arcsec \
        -weight briggs 0.5 \
        -auto-mask 3 \
        -auto-threshold 1 \
        -temp-dir ${tempdir} \
        -join-channels \
        -channels-out $nch \
        -fit-spectral-pol 4 \
        "${tempdir}/${msbase}" | tee "wsclean_${mshead}_${datacolumn}_${target}.log"
    /data/curtin_mwaeor/src/src-workloads/tasks/karabo/scripts/create_cube_with_beam.py ${tempdir}/${mshead}_${datacolumn}_${target} XX-image $nch
    /data/curtin_mwaeor/src/src-workloads/tasks/karabo/scripts/create_cube_with_beam.py ${tempdir}/${mshead}_${datacolumn}_${target} YY-image $nch
    mkdir -p ./img
    mv ${tempdir}/${mshead}_${datacolumn}_${target}*cube.fits ./img
}


declare -A targets
# targets[3C298]="14h19m08.52539063s 6d28m32.74452209s"
targets[HydA]="9h18m04.5252602876s -12d02m04.7299768852s"

MS_DATA=prep/1274010302_2s_40kHz.ms
MS_MODEL=cal/1274010302_2s_40kHz_dp3nosmear_model.ms
for target in "${!targets[@]}"; do
    echo $target
    chg=${targets[$target]};
    run_wsclean "$chg" "${MS_MODEL}" "DATA" "$target"
    run_wsclean "$chg" "${MS_MODEL}" "MODEL_DATA" "$target"
    run_wsclean "$chg" "${MS_DATA}" "CORRECTED_DATA" "$target"
done

exit