From 18513e77fd61ff37c0fa8837582c68ad8e161a61 Mon Sep 17 00:00:00 2001 From: d3v-null Date: Sun, 18 Feb 2024 11:24:22 +0000 Subject: [PATCH] [WIP] recovered from dockerfile after Gacrux crash --- Dockerfile | 42 ++++ src/cli/common/beam/mod.rs | 22 +- src/cli/common/input_vis/mod.rs | 31 ++- src/cli/peel/mod.rs | 10 +- src/gpu/peel.cu | 12 +- src/hyperdrive-peel.code-workspace | 20 ++ src/io/read/uvfits/mod.rs | 5 + src/math/mod.rs | 7 + src/params/input_vis.rs | 22 ++ src/params/peel/mod.rs | 318 ++++++++++++++++++++++++--- src/solutions/mod.rs | 2 +- src/srclist/types/source_list/mod.rs | 25 +++ 12 files changed, 463 insertions(+), 53 deletions(-) create mode 100644 Dockerfile create mode 100644 src/hyperdrive-peel.code-workspace diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 00000000..a2934e2f --- /dev/null +++ b/Dockerfile @@ -0,0 +1,42 @@ +ARG NVIDIA_VERSION=11.4.3 + +FROM nvidia/cuda:${NVIDIA_VERSION}-devel-ubuntu20.04 + +ENV DEBIAN_FRONTEND="noninteractive" +RUN apt-get update -y && \ + apt-get -y install \ + tzdata \ + build-essential \ + pkg-config \ + cmake \ + curl \ + git \ + lcov \ + fontconfig \ + libfreetype-dev \ + libexpat1-dev \ + libcfitsio-dev \ + libhdf5-dev \ + clang \ + libfontconfig-dev \ + && apt-get clean all \ + && rm -rf /var/lib/apt/lists/* + +ARG RUST_VERSION=1.72 +ARG TARGET_CPU=x86-64 +# example: 70 +ARG CUDA_COMPUTE + +# Get Rust +RUN mkdir -m755 /opt/rust /opt/cargo +ENV RUSTUP_HOME=/opt/rust CARGO_HOME=/opt/cargo PATH=/opt/cargo/bin:$PATH +# set minimal rust version here to use a newer stable version +RUN curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal --default-toolchain=$RUST_VERSION + +ADD . /mwa_hyperdrive +WORKDIR /mwa_hyperdrive +ENV CXX=/usr/bin/g++ +ENV CARGO_BUILD_RUSTFLAGS="-C target-cpu=${TARGET_CPU}" +RUN [ -z "$CUDA_COMPUTE" ] || export HYPERDRIVE_CUDA_COMPUTE=${CUDA_COMPUTE}; \ + cargo install --path . --no-default-features --features=cuda,plotting --locked \ + && cargo clean \ No newline at end of file diff --git a/src/cli/common/beam/mod.rs b/src/cli/common/beam/mod.rs index df0b6df3..d3ba1054 100644 --- a/src/cli/common/beam/mod.rs +++ b/src/cli/common/beam/mod.rs @@ -118,8 +118,15 @@ impl BeamArgs { trace!("Attempting to use delays:"); match &dipole_delays { Delays::Full(d) => { - for row in d.outer_iter() { - trace!("{row}"); + let mut last_row = None; + for (i, row) in d.outer_iter().enumerate() { + if let Some(last_row) = last_row { + if row == last_row { + continue + } + } + trace!("{i:03} {row}"); + last_row = Some(row); } } Delays::Partial(d) => trace!("{d:?}"), @@ -178,8 +185,15 @@ impl BeamArgs { }; if let Some(dipole_gains) = dipole_gains.as_ref() { trace!("Attempting to use dipole gains:"); - for row in dipole_gains.outer_iter() { - trace!("{row}"); + let mut last_row = None; + for (i, row) in dipole_gains.outer_iter().enumerate() { + if let Some(last_row) = last_row { + if row == last_row { + continue + } + } + trace!("{i:03} {row}"); + last_row = Some(row); } // Currently, the only way to have dipole gains other than diff --git a/src/cli/common/input_vis/mod.rs b/src/cli/common/input_vis/mod.rs index 4c83454e..f7e3d1c2 100644 --- a/src/cli/common/input_vis/mod.rs +++ b/src/cli/common/input_vis/mod.rs @@ -832,8 +832,9 @@ impl InputVisArgs { r } }; + let first_obs_ts = obs_context.timestamps.first(); time_printer - .push_line(format!("First obs timestamp: {}", obs_context.timestamps.first()).into()); + .push_line(format!("First obs timestamp: {} | {:.2}", first_obs_ts, first_obs_ts.to_gpst_seconds()).into()); time_printer.push_block(vec![ format!( "Available timesteps: {}", @@ -852,26 +853,34 @@ impl InputVisArgs { ) .into()]; match timesteps_to_use.as_slice() { - [t] => block.push( - format!( - "Only timestamp (GPS): {:.2}", - obs_context.timestamps[*t].to_gpst_seconds() + [t] => { + let only_use_ts = obs_context.timestamps[*t]; + block.push( + format!( + "Only timestamp: {} ({:+.2}s)", + only_use_ts, + only_use_ts.to_gpst_seconds() - first_obs_ts.to_gpst_seconds(), + ) + .into(), ) - .into(), - ), + }, [f, .., l] => { + let first_use_ts = obs_context.timestamps[*f]; block.push( format!( - "First timestamp (GPS): {:.2}", - obs_context.timestamps[*f].to_gpst_seconds() + "First timestamp: {} ({:+.2}s)", + first_use_ts, + first_use_ts.to_gpst_seconds() - first_obs_ts.to_gpst_seconds(), ) .into(), ); + let last_use_ts = obs_context.timestamps[*l]; block.push( format!( - "Last timestamp (GPS): {:.2}", - obs_context.timestamps[*l].to_gpst_seconds() + "Last timestamp : {} ({:+.2}s)", + last_use_ts, + last_use_ts.to_gpst_seconds() - first_obs_ts.to_gpst_seconds(), ) .into(), ); diff --git a/src/cli/peel/mod.rs b/src/cli/peel/mod.rs index 7a3708be..28a998b8 100644 --- a/src/cli/peel/mod.rs +++ b/src/cli/peel/mod.rs @@ -28,6 +28,7 @@ use crate::{ params::{ModellingParams, PeelParams}, unit_parsing::{parse_wavelength, WavelengthUnit, WAVELENGTH_FORMATS}, HyperdriveError, + math::div_ceil, }; const DEFAULT_OUTPUT_PEEL_FILENAME: &str = "hyperdrive_peeled.uvfits"; @@ -481,8 +482,8 @@ impl PeelArgs { let n_low_freqs = low_res_spw.get_all_freqs().len(); let n_input_freqs = input_vis_params.spw.get_all_freqs().len(); assert_eq!( - n_low_freqs * low_res_freq_average_factor.get(), - n_input_freqs, + n_low_freqs, + div_ceil(n_input_freqs, low_res_freq_average_factor.get()), "low chans (flagged+unflagged) {} * low_res_freq_average_factor {} != input chans (flagged+unflagged) {}. also iono_freq_average_factor {}", n_low_freqs, low_res_freq_average_factor.get(), @@ -687,12 +688,15 @@ impl PeelArgs { block.push( format!( "- low: {} kHz (averaging {}x)", - low_res_spw.freq_res / 1e3 * low_res_freq_average_factor.get() as f64, + low_res_spw.freq_res / 1e3, low_res_freq_average_factor ) .into(), ); } + block.push( + format!("- {num_passes} passes of {num_loops} loops at {convergence:.2} convergence").into() + ); peel_printer.push_block(block); } diff --git a/src/gpu/peel.cu b/src/gpu/peel.cu index 247d1c23..1fd191f9 100644 --- a/src/gpu/peel.cu +++ b/src/gpu/peel.cu @@ -467,11 +467,12 @@ __global__ void reduce_freqs( // warning: gain is actually needed for model real to line up with vis real. iono_consts->gain *= 1. + convergence * (s_vm / s_mm - 1); // printf( - // "iono: a %+6.4e b %+6.4e a_uu %+6.4e a_uv %+6.4e a_vv %+6.4e A_u %+6.4e A_v %+6.4e denom %+6.4e s_vm %+6.4e s_mm %+6.4e s_vm/s_mm %+6.4e \n", + // "iono: a %+6.4e b %+6.4e a_uu %+6.4e a_uv %+6.4e a_vv %+6.4e A_u %+6.4e A_v %+6.4e denom %+6.4e s_vm %+6.4e s_mm %+6.4e s_vm/s_mm %+6.4e nf %d \n", // iono_consts->alpha, iono_consts->beta, // a_uu, a_uv, a_vv, A_u, A_v, denom, // // iono_consts->gain, - // s_vm, s_mm, s_vm / s_mm + // s_vm, s_mm, s_vm / s_mm, + // num_freqs // ); } } @@ -820,9 +821,10 @@ extern "C" const char *iono_loop(const JonesF32 *d_vis_residual, const float *d_ // Sane? gpuMemcpy(iono_consts, d_iono_consts, sizeof(IonoConsts), gpuMemcpyDeviceToHost); if (iono_consts->gain < 0.0) { - iono_consts->alpha = 0.0; - iono_consts->beta = 0.0; - iono_consts->gain = 1.0; + // let's handle this in cpu world. This gets overwritten by the gpuMemcpy at the end of this function anyway + // iono_consts->alpha = 0.0; + // iono_consts->beta = 0.0; + // iono_consts->gain = 1.0; break; } // printf("iter %d: %.4e %.4e %.4e\n", iteration, iono_consts->alpha, iono_consts->beta, iono_consts->gain); diff --git a/src/hyperdrive-peel.code-workspace b/src/hyperdrive-peel.code-workspace new file mode 100644 index 00000000..16b3338a --- /dev/null +++ b/src/hyperdrive-peel.code-workspace @@ -0,0 +1,20 @@ +{ + "folders": [ + { + "path": ".." + }, + { + "path": "../../../../../data/dev/1226233968" + }, + { + "path": "../../../../../data/dev/1090701368" + }, + { + "path": "../../../../../data/dev/1094090416" + }, + { + "path": "../../../../../data/dev/1096811152" + } + ], + "settings": {} +} \ No newline at end of file diff --git a/src/io/read/uvfits/mod.rs b/src/io/read/uvfits/mod.rs index 1e197f88..d1a49523 100644 --- a/src/io/read/uvfits/mod.rs +++ b/src/io/read/uvfits/mod.rs @@ -1096,6 +1096,11 @@ impl VisRead for UvfitsReader { weights_fb, tile_baseline_flags, }); + debug!( + "[uvfits read_crosses] npols {:?} fppol {:?}", + self.metadata.pols.num_pols(), + self.metadata.num_floats_per_pol + ); match ( self.metadata.pols.num_pols(), self.metadata.num_floats_per_pol, diff --git a/src/math/mod.rs b/src/math/mod.rs index bc574245..90568d07 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -42,7 +42,14 @@ pub(crate) fn average_epoch>(es: I) -> Epoch { Epoch::from_gpst_seconds(average).round(10.milliseconds()) } +// TODO (dev): a.div_ceil(b) would be better, but it's nightly: +// https://doc.rust-lang.org/std/primitive.i32.html#method.div_ceil +pub(crate) fn div_ceil(a: usize, b: usize) -> usize { + (a + b - 1) / b +} + /// Information on flagged tiles, baselines and maps to and from array indices. +#[derive(Debug)] pub struct TileBaselineFlags { /// Map between a pair of tile numbers and its unflagged *cross-correlation* /// baseline index. e.g. If tiles 0 and 2 are flagged, (1, 3) maps to 0 diff --git a/src/params/input_vis.rs b/src/params/input_vis.rs index ec823806..52e507eb 100644 --- a/src/params/input_vis.rs +++ b/src/params/input_vis.rs @@ -65,6 +65,26 @@ pub(crate) struct InputVisParams { pub(crate) dut1: Duration, } +// derive debug for InputVisParams +impl std::fmt::Debug for InputVisParams { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("InputVisParams") + .field("vis_reader.get_input_data_type()", &self.vis_reader.get_input_data_type()) + .field("solutions.is_some()", &self.solutions.is_some()) + .field("timeblocks.timestamps", &self.timeblocks.iter().map(|t| t.timestamps.clone()).collect::>()) + .field("time_res.to_seconds()", &self.time_res.to_seconds()) + .field("spw.chanblocks.len()", &self.spw.chanblocks.len()) + .field("spw.chans_per_chanblock", &self.spw.chans_per_chanblock) + .field("spw.first_freq", &self.spw.first_freq) + .field("spw.flagged_chan_indices", &self.spw.flagged_chan_indices) + .field("tile_baseline_flags.flagged_tiles", &self.tile_baseline_flags.flagged_tiles) + .field("using_autos", &self.using_autos) + .field("ignore_weights", &self.ignore_weights) + .field("dut1.to_seconds()", &self.dut1.to_seconds()) + .finish() + } +} + impl InputVisParams { pub(crate) fn get_obs_context(&self) -> &ObsContext { self.vis_reader.get_obs_context() @@ -103,6 +123,8 @@ impl InputVisParams { assert_eq!(auto_weights_fb.dim(), avg_auto_vis_shape); } + debug!("[read_timeblock] {:?}", self); + let averaging = timeblock.timestamps.len() > 1 || self.spw.chans_per_chanblock.get() > 1; if averaging { diff --git a/src/params/peel/mod.rs b/src/params/peel/mod.rs index d2d81330..45d95b56 100644 --- a/src/params/peel/mod.rs +++ b/src/params/peel/mod.rs @@ -13,6 +13,7 @@ use std::{ ops::{Div, Neg, Sub}, path::PathBuf, thread::{self, ScopedJoinHandle}, + cmp::Ordering, }; use crossbeam_channel::{bounded, unbounded, Receiver, Sender}; @@ -21,7 +22,7 @@ use hifitime::{Duration, Epoch}; use indexmap::IndexMap; use indicatif::{MultiProgress, ProgressBar, ProgressDrawTarget, ProgressStyle}; use itertools::{izip, Itertools}; -use log::{debug, info, trace}; +use log::{debug, info, trace, warn}; use marlu::{ constants::VEL_C, pos::xyz::xyzs_to_cross_uvws, @@ -48,6 +49,7 @@ use crate::{ model::{ModelDevice, ModelError, SkyModeller, SkyModellerCpu}, srclist::SourceList, Chanblock, TileBaselineFlags, MODEL_DEVICE, PROGRESS_BARS, + math::div_ceil, }; #[cfg(any(feature = "cuda", feature = "hip"))] use crate::{ @@ -83,6 +85,48 @@ pub(crate) struct SourceIonoConsts { // pub(crate) centroid_timestamps: Vec, } +#[derive(Debug, PartialEq)] //, Serialize, Deserialize)] +pub(crate) struct BadSource { + // timeblock: Timeblock, + pub(crate) gpstime: f64, + pub(crate) pass: usize, + // pub(crate) gsttime: Epoch, + // pub(crate) times: Vec, + // source: Source, + pub(crate) i_source: usize, + pub(crate) source_name: String, + // pub(crate) weighted_catalogue_pos_j2000: RADec, + // iono_consts: IonoConsts, + pub(crate) alpha: f64, + pub(crate) beta: f64, + pub(crate) gain: f64, + // pub(crate) alphas: Vec, + // pub(crate) betas: Vec, + // pub(crate) gains: Vec, + + pub(crate) residuals_i: Vec>, + pub(crate) residuals_q: Vec>, + pub(crate) residuals_u: Vec>, + pub(crate) residuals_v: Vec>, +} + +// custom sorting implementations +impl PartialOrd for BadSource { + fn partial_cmp(&self, other: &Self) -> Option { + match self.source_name.partial_cmp(&other.source_name) { + // Some(Ordering::Equal) => match self.gpstime.partial_cmp(&other.gpstime) { + Some(Ordering::Equal) => match self.pass.partial_cmp(&other.pass) { + Some(Ordering::Less) => Some(Ordering::Greater), + Some(Ordering::Greater) => Some(Ordering::Less), + other => other, + }, + // other => other, + // }, + other => other, + } + } +} + pub(crate) struct PeelParams { pub(crate) input_vis_params: InputVisParams, pub(crate) output_vis_params: Option, @@ -179,8 +223,7 @@ impl PeelParams { }; assert!( - input_vis_params.timeblocks.len() - == iono_time_average_factor.get() * iono_timeblocks.len(), + div_ceil(input_vis_params.timeblocks.len(), iono_time_average_factor.get()) == iono_timeblocks.len(), "num_read_times {} != num_iono_times {} * iono_time_average_factor {}", input_vis_params.timeblocks.len(), iono_timeblocks.len(), @@ -188,8 +231,8 @@ impl PeelParams { ); let error = AtomicCell::new(false); - let (tx_data, rx_data) = bounded(1); - let (tx_residual, rx_residual) = bounded(1); + let (tx_data, rx_data) = bounded(2); + let (tx_residual, rx_residual) = bounded(2); let (tx_full_residual, rx_full_residual) = bounded(iono_time_average_factor.get()); let (tx_write, rx_write) = bounded(2); let (tx_iono_consts, rx_iono_consts) = unbounded(); @@ -513,12 +556,6 @@ fn get_weights_rts(tile_uvs: ArrayView2, lambdas_m: &[f64], short_sigma: f64 weights } -// TODO (dev): a.div_ceil(b) would be better, but it's nightly: -// https://doc.rust-lang.org/std/primitive.i32.html#method.div_ceil -fn div_ceil(a: usize, b: usize) -> usize { - (a + b - 1) / b -} - /// Like `vis_weight_average_tfb`, but for when we don't need to keep the low-res weights /// Average "high-res" vis and weights to "low-res" vis (no low-res weights) /// arguments are all 3D arrays with axes (time, freq, baseline). @@ -1428,6 +1465,10 @@ fn peel_gpu( no_precession: bool, multi_progress_bar: &MultiProgress, ) -> Result<(), PeelError> { + use std::collections::{HashMap, HashSet}; + + use crate::srclist::{ComponentType, FluxDensityType, FluxDensity}; + let timestamps = &timeblock.timestamps; let num_timesteps = vis_residual_tfb.len_of(Axis(0)); @@ -1490,6 +1531,8 @@ fn peel_gpu( ); peel_progress.tick(); + macro_rules! pb_warn { ($($arg:tt)+) => (multi_progress_bar.suspend(|| warn!($($arg)+))) } + macro_rules! pb_info { ($($arg:tt)+) => (multi_progress_bar.suspend(|| info!($($arg)+))) } macro_rules! pb_debug { ($($arg:tt)+) => (multi_progress_bar.suspend(|| debug!($($arg)+))) } macro_rules! pb_trace { ($($arg:tt)+) => (multi_progress_bar.suspend(|| trace!($($arg)+))) } @@ -1653,8 +1696,10 @@ fn peel_gpu( .try_into() .expect("smaller than i32::MAX"); + let mut bad_sources = Vec::::new(); + unsafe { - let mut gpu_uvws = Array2::default((num_timesteps, num_cross_baselines)); + let mut gpu_uvws: ArrayBase, Dim<[usize; 2]>> = Array2::default((num_timesteps, num_cross_baselines)); gpu_uvws .outer_iter_mut() .zip(tile_xyzs_high_res.outer_iter()) @@ -1774,6 +1819,11 @@ fn peel_gpu( pass + 1 )); + let mut pass_issues = 0; + let mut pass_alpha_mag = 0.; + let mut pass_bega_mag = 0.; + let mut pass_gain_mag = 0.; + // this needs to be inside the pass loop, because d_uvws_from gets swapped with d_uvws_to gpu_kernel_call!( gpu::xyzs_to_uvws, @@ -1797,9 +1847,9 @@ fn peel_gpu( .enumerate() { let start = std::time::Instant::now(); - pb_debug!( - "peel loop {pass}: {source_name} at {source_pos} (has iono {iono_consts:?})" - ); + // pb_debug!( + // "peel loop {pass}: {source_name} at {source_pos} (has iono {iono_consts:?})" + // ); // let old_iono_consts = *iono_consts; let old_iono_consts = IonoConsts { @@ -1944,7 +1994,7 @@ fn peel_gpu( // we have low res residual and model, now print what iono fit will see. // { // use marlu::math::cross_correlation_baseline_to_tiles; - // let vis_residual_low_res_fb = d_low_res_resid_fb.copy_from_device_new()?; + // let vis_residual_low_res_fb = d_low_res_vis_fb.copy_from_device_new()?; // dbg!(vis_residual_low_res_fb.len(), num_low_res_chans, num_cross_baselines); // let vis_residual_low_res_fb = Array2::from_shape_vec( // (num_low_res_chans, num_cross_baselines), @@ -1992,7 +2042,7 @@ fn peel_gpu( // let model_i = model[0] + model[3]; // let u = u as GpuFloat; // let v = v as GpuFloat; - // println!("uv {ant1:3} {ant2:3} ({u:+9.3}, {v:+9.3}) l{lambda:+7.5} wt{weight:3.1} | RI {:+11.7} @{:+5.3}pi | MI {:+11.7} @{:+5.3}pi", residual_i.norm(), residual_i.arg(), model_i.norm(), model_i.arg()); + // println!("uv {ant1:3} {ant2:3} ({u:+9.3}, {v:+9.3}) l{lambda:+7.5} wt{weight:+3.1} | RI {:+11.7} @{:+5.3}pi | MI {:+11.7} @{:+5.3}pi", residual_i.norm(), residual_i.arg(), model_i.norm(), model_i.arg()); // } // } // } @@ -2002,6 +2052,18 @@ fn peel_gpu( beta: 0.0, gain: 1.0, }; + // get size of device ptr + let lrblch = (num_cross_baselines_i32 * num_low_res_chans_i32) as f64; + pb_trace!("before iono_loop nt{:?} nxbl{:?} nlrch{:?} = lrxblch{:?}; lrvfb{:?} lrwfb{:?} lrmfb{:?} lrmrfb{:?}", + num_tiles_i32, + num_cross_baselines_i32, + num_low_res_chans_i32, + lrblch, + d_low_res_vis_fb.get_size() as f64 / lrblch, + d_low_res_weights_fb.get_size() as f64 / lrblch, + d_low_res_model_fb.get_size() as f64 / lrblch, + d_low_res_model_rotated.get_size() as f64 / lrblch, + ); gpu_kernel_call!( gpu::iono_loop, d_low_res_vis_fb.get().cast(), @@ -2024,6 +2086,76 @@ fn peel_gpu( iono_consts.alpha = old_iono_consts.alpha + gpu_iono_consts.alpha; iono_consts.beta = old_iono_consts.beta + gpu_iono_consts.beta; iono_consts.gain = old_iono_consts.gain * gpu_iono_consts.gain; + + let issues = format!( + "{}{}{}", + if iono_consts.alpha.abs() > 1e-3 { + if iono_consts.alpha > 0.0 { "A" } else { "a" } + } else { + "" + }, + if iono_consts.beta.abs() > 1e-3 { + if iono_consts.beta > 0.0 { "B" } else { "b" } + } else { + "" + }, + if iono_consts.gain < 0.0 { + "g" + } else if iono_consts.gain > 1.5 { + "G" + } else { + "" + }, + ); + let message = format!( + // "t{:03} p{pass} s{i_source:6}|{source_name:16} @ ra {:+7.2} d {:+7.2} | a {:+7.2e} b {:+7.2e} g {:+3.2} | da {:+8.2e} db {:+8.2e} dg {:+3.2} | {}", + "t{:3} pass {:2} s{i_source:6}|{source_name:16} @ ra {:+7.2} d {:+7.2} | a {:+8.6} b {:+8.6} g {:+3.2} | da {:+8.6} db {:+8.6} dg {:+3.2} | {}", + timeblock.index, + pass+1, + (source_pos.ra.to_degrees() + 180.) % 360. - 180., + source_pos.dec.to_degrees(), + iono_consts.alpha, + iono_consts.beta, + iono_consts.gain, + iono_consts.alpha - old_iono_consts.alpha, + iono_consts.beta - old_iono_consts.beta, + iono_consts.gain - old_iono_consts.gain, + issues, + ); + if issues.is_empty() { + pb_debug!("[peel_gpu] {}", message); + pass_alpha_mag += (iono_consts.alpha - old_iono_consts.alpha).abs(); + pass_bega_mag += (iono_consts.beta - old_iono_consts.beta).abs(); + pass_gain_mag += iono_consts.gain - old_iono_consts.gain; + } else { + pb_debug!( + "[peel_gpu] {} (reverting to a {:+8.6} b {:+8.6} g {:+3.2})", + message, + old_iono_consts.alpha, + old_iono_consts.beta, + old_iono_consts.gain, + ); + bad_sources.push(BadSource { + gpstime: timeblock.median.to_gpst_seconds(), + pass, + i_source, + source_name: source_name.to_string(), + // weighted_catalogue_pos_j2000: source_pos, + alpha: iono_consts.alpha, + beta: iono_consts.beta, + gain: iono_consts.gain, + residuals_i: Vec::default(), // todo!(), + residuals_q: Vec::default(), // todo!(), + residuals_u: Vec::default(), // todo!(), + residuals_v: Vec::default(), // todo!(), + }); + + iono_consts.alpha = old_iono_consts.alpha; + iono_consts.beta = old_iono_consts.beta; + iono_consts.gain = old_iono_consts.gain; + pass_issues += 1; + } + let gpu_iono_consts = gpu::IonoConsts { alpha: iono_consts.alpha, beta: iono_consts.beta, @@ -2097,12 +2229,32 @@ fn peel_gpu( // The new phase centre becomes the old one. std::mem::swap(&mut d_uvws_from, &mut d_uvws_to); - pb_debug!( - "peel loop finished: {source_name} at {source_pos} (has iono {iono_consts:?})" - ); peel_progress.inc(1); } + let num_good_sources = (num_sources_to_iono_subtract - pass_issues) as f64; + if num_good_sources > 0. { + let msg = format!( + "t{:3} pass {:2} ma {:+7.2e} mb {:+7.2e} mg {:+7.2e}", + timeblock.index, + pass + 1, + pass_alpha_mag / num_good_sources, + pass_bega_mag / num_good_sources, + pass_gain_mag / num_good_sources + ); + if pass_issues > 0 { + pb_warn!("[peel_gpu] {} ({} issues)", msg, pass_issues); + } else { + pb_info!("[peel_gpu] {} (no issues)", msg); + } + } else { + pb_warn!( + "[peel_gpu] t{:03} pass {:2} all sources had issues", + timeblock.index, + pass + 1 + ); + } + // Rotate back to the phase centre. gpu_kernel_call!( gpu::xyzs_to_uvws, @@ -2134,6 +2286,100 @@ fn peel_gpu( d_high_res_vis_tfb.copy_from_device(vis_residual_tfb.as_slice_mut().unwrap())?; } + let mut pass_counts = HashMap::::new(); + for bad_source in bad_sources.iter() { + *pass_counts.entry(bad_source.source_name.clone()).or_default() += 1; + } + let mut printed = HashSet::::new(); + bad_sources.sort_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal)); + for bad_source in bad_sources.into_iter() { + let BadSource { + gpstime, + // pass, + i_source, + source_name, + alpha, + beta, + gain, + .. + } = bad_source; + // if source_name is in printed + if printed.contains(&source_name) { + continue; + } else { + printed.insert(source_name.clone()); + } + let passes = pass_counts[&source_name]; + + let pos = source_weighted_positions[i_source]; + let RADec { ra, dec } = pos; + pb_warn!( + "[peel_gpu] Bad source: {:.2} p{:2} {} at radec({:+7.2}, {:+7.2}) iono({:+8.6},{:+8.6},{:+3.2})", + gpstime, + passes, + source_name, + (ra + 180.) % 360. - 180., + dec, + alpha, + beta, + gain + ); + let matches = source_list.search_asec(pos, 300.); + for (sep, src_name, idx, comp) in matches { + let compstr = match comp.comp_type { + ComponentType::Gaussian { maj, min, pa } => format!( + "G {:6.1}as {:6.1}as {:+6.1}d", + maj.to_degrees() * 3600., + min.to_degrees() * 3600., + pa.to_degrees() + ), + ComponentType::Point => format!( + "P {:8} {:8} {:7}", + "", "", "" + ), + _ => format!( + "? {:8} {:8} {:7}", + "", "", "" + ), + }; + let fluxstr = match comp.flux_type { + FluxDensityType::CurvedPowerLaw { si, fd: FluxDensity { freq, i, .. }, q } => format!( + "cpl S={:+6.2}(νn)^{:+5.2} exp[{:+5.2}ln(νn)]; @{:.1}MHz", i, si, q, freq / 1e6 + ), + FluxDensityType::PowerLaw { si, fd: FluxDensity { freq, i, .. } } => format!( + "pl S={:+6.2}(νn)^{:+5.2}; @{:.1}MHz", i, si, freq / 1e6 + ), + FluxDensityType::List(fds) => { + let FluxDensity { i, freq, .. } = fds[0]; + format!("lst S={:+6.2} @{:.1}MHz", i, freq / 1e6) + }, + }; + pb_warn!( + "[peel_gpu] {sep:5.1} {src_name:16} c{idx:2} at radec({:+7.2},{:+7.2}) comp({}) flx ({})", + (comp.radec.ra + 180.) % 360. - 180., + comp.radec.dec, + compstr, + fluxstr, + ); + } + // pb_warn!( + // " residuals_i: {:?}", + // bad_source.residuals_i.iter().map(|c| c.norm()).collect_vec() + // ); + // pb_warn!( + // " residuals_q: {:?}", + // bad_source.residuals_q.iter().map(|c| c.norm()).collect_vec() + // ); + // pb_warn!( + // " residuals_u: {:?}", + // bad_source.residuals_u.iter().map(|c| c.norm()).collect_vec() + // ); + // pb_warn!( + // " residuals_v: {:?}", + // bad_source.residuals_v.iter().map(|c| c.norm()).collect_vec() + // ); + } + Ok(()) } @@ -2279,6 +2525,14 @@ fn read_thread( return Ok(()); } + debug!( + "[read] vdfb shp={:?} sum={:?} vwfb shp={:?} sum={:?}", + vis_data_fb.shape(), + vis_data_fb.sum(), + vis_weights_fb.shape(), + vis_weights_fb.sum(), + ); + match tx_data.send((vis_data_fb, vis_weights_fb, timeblock.median)) { Ok(()) => (), // If we can't send the message, it's because the @@ -2418,7 +2672,8 @@ fn subtract_thread( // residuals. vis_data_fb.iter_mut().for_each(|j| *j *= -1.0); - let (lst, xyzs, latitude) = if apply_precession { + // let (lst, xyzs, latitude) = + if apply_precession { let precession_info = precess_time( array_position.longitude_rad, array_position.latitude_rad, @@ -2489,14 +2744,15 @@ fn subtract_thread( #[cfg(any(feature = "cuda", feature = "hip"))] if let Some(GpuStuff { modeller, - d_uvws_from, - d_uvws_to, - d_beam_jones, - d_xyzs, - gpu_xyzs, - d_lmst, - d_lambdas, - d_vis_fb, + // d_uvws_from, + // d_uvws_to, + // d_beam_jones, + // d_xyzs, + // gpu_xyzs, + // d_lmst, + // d_lambdas, + // d_vis_fb, + .. }) = gpu_modeller.as_mut() { // d_vis_fb.overwrite(vis_data_fb.as_slice().expect("is contiguous"))?; @@ -2657,6 +2913,10 @@ fn joiner_thread<'a>( full_weights_fb.assign(&vis_weights_fb); } + if vis_weights_tfb.sum() < 0.0 { + warn!("[joiner] all flagged: timestamps={timestamps:?}") + } + match tx_full_residual.send((vis_residual_tfb, vis_weights_tfb, timeblock)) { Ok(()) => (), Err(_) => return, diff --git a/src/solutions/mod.rs b/src/solutions/mod.rs index f9d1e1f5..a29a1401 100644 --- a/src/solutions/mod.rs +++ b/src/solutions/mod.rs @@ -47,7 +47,7 @@ pub(crate) enum CalSolutionType { Bin, } -#[derive(Default)] +#[derive(Default, Debug)] pub struct CalibrationSolutions { /// The direction-independent calibration solutions. This has dimensions of /// (num_timeblocks, total_num_tiles, total_num_chanblocks). Note that this diff --git a/src/srclist/types/source_list/mod.rs b/src/srclist/types/source_list/mod.rs index 6d67590b..1058990c 100644 --- a/src/srclist/types/source_list/mod.rs +++ b/src/srclist/types/source_list/mod.rs @@ -13,6 +13,8 @@ use std::ops::{Deref, DerefMut}; use indexmap::IndexMap; use serde::{Deserialize, Serialize}; +use marlu::RADec; + use super::*; /// A [`IndexMap`] of source names for keys and [`Source`] structs for values. @@ -83,6 +85,29 @@ impl SourceList { .collect(); SourceList(sl) } + + /// Find components within a given angular distance [radians] of a position + /// Return a vector of tuples containing separation [radians], source name, component index, and component. + pub fn search(&self, target: RADec, radians: f64) -> Vec<(f64, String, usize, SourceComponent)> { + let mut matches = Vec::new(); + for (name, src) in self.iter() { + for (i, comp) in src.components.iter().enumerate() { + let sep = comp.radec.separation(target); + if sep < radians { + matches.push((sep, name.clone(), i, comp.clone())); + } + } + } + matches + } + + /// Find components within a given angular distance [arcseconds] of a position + /// Return a vector of tuples containing separation [radians], source name, component index, and component. + pub fn search_asec(&self, target: RADec, arcseconds: f64) -> Vec<(f64, String, usize, SourceComponent)> { + self.search(target, arcseconds.to_radians() / 3600.0).into_iter().map(|(sep, name, i, comp)| + (sep.to_degrees() * 3600.0, name, i, comp) + ).collect() + } } impl From> for SourceList {