Parametric Model Document Formulas¶
This notebook attempts to generate the equations as shown in the parametric model document. Useful for checking whether the parametric model as implemented actually matches the document.
In [1]:
import sys
from IPython.display import display, Math, Latex, HTML, Markdown
from sympy import latex, simplify, Lambda, Symbol, Function
sys.path+=['..']
from sdp_par_model import evaluate
from sdp_par_model import reports as iapi
from sdp_par_model.config import PipelineConfig
from sdp_par_model.parameters.definitions import *
from sdp_par_model.parameters.container import BLDep
In [2]:
# Telescope and band should not matter
cfg = PipelineConfig(telescope=Telescopes.SKA1_Mid, band=Bands.Mid1, pipeline=Pipelines.ICAL)
cfg_rcal = PipelineConfig(telescope=Telescopes.SKA1_Mid, band=Bands.Mid1, pipeline=Pipelines.RCAL)
cfg_ingest = PipelineConfig(telescope=Telescopes.SKA1_Mid, band=Bands.Mid1, pipeline=Pipelines.Ingest)
cfg_specfit = PipelineConfig(telescope=Telescopes.SKA1_Mid, band=Bands.Mid1, pipeline=Pipelines.DPrepA_Image)
tp = cfg.calc_tel_params(symbolify='product')
tp2 = cfg.calc_tel_params(symbolify='helper')
tp3 = cfg.calc_tel_params(symbolify='all')
tp_rcal = cfg_rcal.calc_tel_params(symbolify='product')
tp_ingest = cfg_ingest.calc_tel_params(symbolify='product')
tp_ingest2 = cfg_ingest.calc_tel_params(symbolify='helper')
tp_specfit = cfg_specfit.calc_tel_params(symbolify='product')
def show_sympy(*exprs):
l = "=".join([latex(e) for e in exprs])
display(Math("\\begin{aligned}" + l + "\\end{aligned}"))
In [3]:
b = Symbol("b")
tp_ingest.Rviscorr = tp_ingest.Nbeam * tp_ingest.Npp * tp_ingest.Rvis_ingest
tp_ingest2.Rviscorr = tp_ingest2.Nbeam * tp_ingest2.Npp * tp_ingest2.Rvis_ingest
tp.RvisfullFoV = tp_ingest2.Nbeam * tp_ingest2.Npp * tp.Rvis_predict
# Definition of Rvis consistent with most uses in the document, and a helper to introduce it into formulas
RvisfullFoV = tp.Nbeam * tp.Nbl * tp.Npp * tp.Nf_vis / tp.Tint_used
def rvisfullfov(expr):
return expr / RvisfullFoV(b) * Symbol("R_vis,fullFoV")
Imaging¶
In [4]:
# Ingest rate over all polarisations and beams.
show_sympy(Symbol("R_vis,corr"), tp_ingest.Rviscorr, tp_ingest2.Rviscorr)
$\displaystyle \begin{aligned}R_{vis,corr}=N_{beam} N_{pp} R_{vis,ingest}=\frac{N_{beam} N_{f,max} N_{pp} \left(N_{a} + N_{bl}\right)}{T_{int,used}}\end{aligned}$
In [5]:
# Pixel production rate (after reprojection), assumed to be summed over the entire observation
show_sympy(Symbol("R_pix"), tp.products[Products.Reprojection]['Rout'] * tp.Tsnap / tp.Tobs)
$\displaystyle \begin{aligned}R_{pix}=\frac{M_{px} N_{beam} N_{f,proj,backward} N_{facet}^{2} N_{majortotal} N_{pix,linear}^{2} N_{pp}}{T_{obs}}\end{aligned}$
Receive¶
In [6]:
show_sympy(Symbol("C_Receive"), tp_ingest.products[Products.Receive]['Rflop'])
$\displaystyle \begin{aligned}C_{Receive}=N_{beam} N_{f,min} \left(\frac{1000 N_{a}}{T_{int,used}} + \frac{4 N_{pp} R_{vis,ingest}}{N_{f,min}}\right)\end{aligned}$
Flagging¶
In [7]:
show_sympy(Symbol("N_flop/vis"), tp_ingest.products[Products.Flag]['Rflop'] / tp_ingest.Rviscorr)
show_sympy(Symbol("C_Flag"), tp_ingest.products[Products.Flag]['Rflop'])
$\displaystyle \begin{aligned}N_{flop/vis}=278\end{aligned}$
$\displaystyle \begin{aligned}C_{Flag}=278 N_{beam} N_{pp} R_{vis,ingest}\end{aligned}$
Demix¶
In [8]:
show_sympy(Symbol("N_flop/vis"), tp_ingest.products[Products.Demix]['Rflop'] / tp_ingest.Rviscorr)
show_sympy(Symbol("C_Demix"), tp_ingest.products[Products.Demix]['Rflop'])
$\displaystyle \begin{aligned}N_{flop/vis}=\frac{N_{f,min} \left(154 N_{Ateam} + 84 + \frac{N_{f,min} T_{int,used} \left(N_{Ateam}^{2} + N_{Ateam} \left(24 N_{solve} + 33\right) + 64\right)}{N_{f,max} T_{ion}}\right)}{N_{pp}}\end{aligned}$
$\displaystyle \begin{aligned}C_{Demix}=N_{beam} N_{f,min} R_{vis,ingest} \left(154 N_{Ateam} + 84 + \frac{N_{f,min} T_{int,used} \left(N_{Ateam}^{2} + N_{Ateam} \left(24 N_{solve} + 33\right) + 64\right)}{N_{f,max} T_{ion}}\right)\end{aligned}$
Averaging¶
In [9]:
show_sympy(Symbol("C_Average"), tp_ingest.products[Products.Average]['Rflop'])
$\displaystyle \begin{aligned}C_{Average}=8 N_{beam} N_{pp} R_{vis,ingest}\end{aligned}$
Predict via Direct Fourier Transform¶
In [10]:
show_sympy(Symbol("N_flop/vis^predict,DFT"), tp.products[Products.DFT]['Rflop'] / RvisfullFoV(b) / tp.Nmajortotal)
show_sympy(Symbol("C_Predict"), rvisfullfov(tp.products[Products.DFT]['Rflop']))
$\displaystyle \begin{aligned}N^{predict,DFT}_{flop/vis}=\frac{T_{int,used} \left(32 N_{a}^{2} N_{source} + 128 N_{a}^{2} + 266 N_{a} N_{source}\right) \sum_{i=1}^{N_{bl}} R_{vis}{\left(B_{max}{\left(i \right)},1 \right)}}{N_{bl}^{2} N_{pp} N_{f,vis}{\left(b \right)}}\end{aligned}$
$\displaystyle \begin{aligned}C_{Predict}=\frac{N_{majortotal} R_{vis,fullFoV} T_{int,used} \left(32 N_{a}^{2} N_{source} + 128 N_{a}^{2} + 266 N_{a} N_{source}\right) \sum_{i=1}^{N_{bl}} R_{vis}{\left(B_{max}{\left(i \right)},1 \right)}}{N_{bl}^{2} N_{pp} N_{f,vis}{\left(b \right)}}\end{aligned}$
Solve¶
In [11]:
show_sympy(Symbol("N_{flop/solution}^{StEFCal}"), tp.NFlop_solver)
show_sympy(Symbol("C_{Solve}^{RCAL}"), tp_rcal.products[Products.Solve]['Rflop'])
show_sympy(Symbol("C_{Solve}^{ICAL}"), tp.products[Products.Solve]['Rflop'])
$\displaystyle \begin{aligned}N_{flop/solution}^{StEFCal}=48 N_{a}^{2} N_{pp} N_{solve}\end{aligned}$
$\displaystyle \begin{aligned}C_{Solve}^{RCAL}=N_{beam} N_{subbands} \left(N_{\mathcal{self}} + 1\right) \left(\frac{48 N_{a}^{2} N_{pp} N_{solve} \max\left(1, \frac{N_{f,out}}{N_{subbands}}\right) \max\left(1, \frac{T_{solve}}{t_{RCAL,G}}\right)}{T_{solve}} + \frac{40 N_{pp} \sum_{i=1}^{N_{bl}} R_{vis}{\left(B_{max}{\left(i \right)},1 \right)}}{N_{subbands}}\right)\end{aligned}$
$\displaystyle \begin{aligned}C_{Solve}^{ICAL}=N_{beam} N_{subbands} \left(N_{\mathcal{self}} + 1\right) \left(\frac{48 N_{a}^{2} N_{pp} N_{solve} \left(N_{Ipatches} \max\left(1, \frac{T_{solve}}{t_{ICAL,I}}\right) + \max\left(1, \frac{N_{B,parameters}}{N_{subbands}}\right) \max\left(1, \frac{T_{solve}}{t_{ICAL,B}}\right) + \max\left(1, \frac{T_{solve}}{t_{ICAL,G}}\right)\right)}{T_{solve}} + \frac{40 N_{pp} \left(N_{Ipatches} + 2\right) \sum_{i=1}^{N_{bl}} R_{vis}{\left(B_{max}{\left(i \right)},1 \right)}}{N_{subbands}}\right)\end{aligned}$
Subtract¶
In [12]:
show_sympy(Symbol("C_Subtract"), tp.products[Products.Subtract_Visibility]['Rflop'])
$\displaystyle \begin{aligned}C_{Subtract}=2 N_{beam} N_{majortotal} N_{pp} \sum_{i=1}^{N_{bl}} R_{vis}{\left(B_{max}{\left(i \right)},1 \right)}\end{aligned}$
Correct¶
In [13]:
show_sympy(Symbol("C_Correct"), tp.products[Products.Correct]['Rflop'])
$\displaystyle \begin{aligned}C_{Correct}=8 N_{Ipatches} N_{beam} N_{majortotal} N_{mm} N_{pp} \sum_{i=1}^{N_{bl}} R_{vis}{\left(B_{max}{\left(i \right)},1 \right)}\end{aligned}$
Phase Rotation¶
In [14]:
show_sympy(Symbol("C_phrot,back"), simplify(tp.products[Products.PhaseRotation]['Rflop']), RvisfullFoV(b))
show_sympy(Symbol("C_phrot,predict"), rvisfullfov(tp.products[Products.PhaseRotationPredict]['Rflop']))
$\displaystyle \begin{aligned}C_{phrot,back}=\begin{cases} 0 & \text{for}\: N_{facet} = 1 \\\frac{28 N_{beam} N_{facet}^{2} N_{majortotal} N_{pp} \left(N_{facet} - 1\right) \sum_{i=1}^{N_{bl}} N_{f,vis}{\left(B_{max}{\left(i \right)} \right)}}{T_{int,used} \left|{N_{facet} - 1}\right|} & \text{otherwise} \end{cases}=\frac{N_{beam} N_{bl} N_{pp} N_{f,vis}{\left(b \right)}}{T_{int,used}}\end{aligned}$
$\displaystyle \begin{aligned}C_{phrot,predict}=\frac{28 N_{facet}^{2} N_{majortotal} R_{vis,fullFoV} \operatorname{sign}{\left(N_{facet} - 1 \right)} \sum_{i=1}^{N_{bl}} N_{f,vis}{\left(B_{max}{\left(i \right)} \right)}}{N_{bl} N_{f,vis}{\left(b \right)}}\end{aligned}$
Grid and Degrid¶
In [15]:
show_sympy(Symbol("C_{grid}^{ICAL}"), tp.products[Products.Grid]['Rflop'])
show_sympy(Symbol("C_{degrid}^{ICAL,A}"), tp.products[Products.Degrid]['Rflop'])
$\displaystyle \begin{aligned}C_{grid}^{ICAL}=8 N_{beam} N_{facet}^{2} N_{majortotal} N_{mm} N_{pp} \sum_{i=1}^{N_{bl}} \frac{N_{f,vis,backward}{\left(B_{max}{\left(i \right)} \right)} N_{kernel2,backward}{\left(B_{max}{\left(i \right)} \right)}}{T_{coal,backward}{\left(B_{max}{\left(i \right)} \right)}}\end{aligned}$
$\displaystyle \begin{aligned}C_{degrid}^{ICAL,A}=8 N_{beam} N_{facet,predict}^{2} N_{majortotal} N_{mm} N_{pp} \sum_{i=1}^{N_{bl}} \frac{N_{f,vis,predict}{\left(B_{max}{\left(i \right)} \right)} N_{kernel2,predict}{\left(B_{max}{\left(i \right)} \right)}}{T_{coal,predict}{\left(B_{max}{\left(i \right)} \right)}}\end{aligned}$
Gridding Kernel Update¶
In [16]:
show_sympy(tp.Nkernel_AW_predict(b), tp2.Nkernel_AW_predict.term)
show_sympy(Function("N_cvff")(b), (tp.Ngcf_used_backward * tp.Nkernel_AW_predict).term)
show_sympy(Function("C_{Kernels,back}"), tp.products[Products.Gridding_Kernel_Update]['Rflop'])
show_sympy(Function("C_{Kernels,predict}"), tp.products[Products.Degridding_Kernel_Update]['Rflop'])
$\displaystyle \begin{aligned}N_{kernel,AW,predict}{\left(b \right)}=\left(N_{aa}^{2} + 1.0 \Theta_{fov,predict}^{2} \left(\Theta_{fov,predict}^{2} \min\left(\Delta_{W stack}, \Delta_{W max}{\left(b \right)}\right)^{2} + \frac{0.636619772367581 \Theta_{fov,predict} \min\left(\Delta_{W stack}, \Delta_{W max}{\left(b \right)}\right)^{1.5}}{\epsilon_{w}}\right)\right)^{0.5}\end{aligned}$
$\displaystyle \begin{aligned}N_{cvff}{\left(b \right)}=N_{gcf,used,backward}{\left(b \right)} N_{kernel,AW,predict}{\left(b \right)}\end{aligned}$
$\displaystyle \begin{aligned}C_{Kernels,back}=\frac{5.0 N_{beam} N_{facet}^{2} N_{majortotal} N_{mm} N_{pp} \sum_{i=1}^{N_{bl}} \frac{N_{f,gcf,backward}{\left(B_{max}{\left(i \right)} \right)} N_{gcf,used,backward}{\left(B_{max}{\left(i \right)} \right)} N_{kernel,AW,backward}^{2}{\left(B_{max}{\left(i \right)} \right)} \log{\left(N_{kernel,AW,backward}^{2}{\left(B_{max}{\left(i \right)} \right)} \right)}}{T_{kernel,backward}{\left(B_{max}{\left(i \right)} \right)}}}{\log{\left(2 \right)}}\end{aligned}$
$\displaystyle \begin{aligned}C_{Kernels,predict}=\frac{5.0 N_{beam} N_{facet,predict}^{2} N_{majortotal} N_{mm} N_{pp} \sum_{i=1}^{N_{bl}} \frac{N_{f,gcf,predict}{\left(B_{max}{\left(i \right)} \right)} N_{gcf,used,predict}{\left(B_{max}{\left(i \right)} \right)} N_{kernel,AW,predict}^{2}{\left(B_{max}{\left(i \right)} \right)} \log{\left(N_{kernel,AW,predict}^{2}{\left(B_{max}{\left(i \right)} \right)} \right)}}{T_{kernel,predict}{\left(B_{max}{\left(i \right)} \right)}}}{\log{\left(2 \right)}}\end{aligned}$
FFT and iFFT¶
In [17]:
show_sympy(tp.Nf_FFT_backward, tp3.Nf_FFT_backward)
show_sympy(Symbol("C_{FFT}^{ICAL}"), tp.products[Products.FFT]['Rflop'])
show_sympy(tp.Nf_FFT_predict, tp3.Nf_FFT_predict)
show_sympy(Symbol("C_{IFFT}^{ICAL,A}"), tp.products[Products.IFFT]['Rflop'])
$\displaystyle \begin{aligned}N_{f,FFT,backward}=N_{f,out}\end{aligned}$
$\displaystyle \begin{aligned}C_{FFT}^{ICAL}=\frac{2.5 N_{beam} N_{f,FFT,backward} N_{facet}^{2} N_{majortotal} N_{pix,linear}^{2} N_{pp} \log{\left(N_{pix,linear}^{2} \right)} \max\left(1, \frac{\Delta_{W max}{\left(B_{max} \right)}}{\Delta_{W stack}}\right)}{T_{snap} \log{\left(2 \right)}}\end{aligned}$
$\displaystyle \begin{aligned}N_{f,FFT,predict}=\min\left(N_{f,max}, N_{f,out} N_{tt}\right)\end{aligned}$
$\displaystyle \begin{aligned}C_{IFFT}^{ICAL,A}=\frac{2.5 N_{beam} N_{f,FFT,predict} N_{facet,predict}^{2} N_{majortotal} N_{pix,linear,predict}^{2} N_{pp} \log{\left(N_{pix,linear,predict}^{2} \right)} \max\left(1, \frac{\Delta_{W max}{\left(B_{max} \right)}}{\Delta_{W stack}}\right)}{T_{snap} \log{\left(2 \right)}}\end{aligned}$
Reprojection¶
In [18]:
show_sympy(Symbol("C_{Reproj}^{ICAL}"), tp.products[Products.Reprojection]['Rflop'])
$\displaystyle \begin{aligned}C_{Reproj}^{ICAL}=\frac{50.0 N_{beam} N_{f,proj,backward} N_{facet}^{2} N_{majortotal} N_{pix,linear}^{2} N_{pp}}{T_{snap}}\end{aligned}$
Image Spectral Fitting¶
In [19]:
show_sympy(Symbol("C_{SpecFit}^{ICAL}"), tp_specfit.products[Products.Image_Spectral_Fitting]['Rflop'])
$\displaystyle \begin{aligned}C_{SpecFit}^{ICAL}=\frac{N_{beam} N_{majortotal} N_{pix,linear,fov,total}^{2} N_{pp} N_{tt} \left(2.0 N_{f,FFT,backward} + 2.0 N_{f,FFT,predict}\right)}{T_{obs}}\end{aligned}$
Subtract Image Component¶
In [20]:
show_sympy(Symbol("C_{grid}^{ICAL}"), tp.products[Products.Identify_Component]['Rflop'])
show_sympy(Symbol("C_{degrid}^{ICAL,A}"), tp.products[Products.Subtract_Image_Component]['Rflop'])
$\displaystyle \begin{aligned}C_{grid}^{ICAL}=\frac{N_{beam} N_{facet}^{2} N_{majortotal} N_{scales} N_{subbands} N_{tt} \left(2 N_{minor} N_{patch}^{2} + 2 N_{pix,linear}^{2}\right)}{T_{obs}}\end{aligned}$
$\displaystyle \begin{aligned}C_{degrid}^{ICAL,A}=\frac{2 N_{beam} N_{majortotal} N_{minor} N_{patch}^{2} N_{pp} N_{scales} N_{subbands} N_{tt}}{T_{obs}}\end{aligned}$
Source find¶
In [21]:
show_sympy(Symbol("C_{Source Find}"), tp.products[Products.Source_Find]['Rflop'])
$\displaystyle \begin{aligned}C_{Source Find}=\frac{600 N_{majortotal} N_{source} N_{source,find,iterations}}{T_{obs}}\end{aligned}$
Time Smearing Limit, tsmear¶
In [22]:
show_sympy(Function("t_{smear}")(b), tp2.Tint_used * tp2.combine_time_samples_facet.term)
$\displaystyle \begin{aligned}t_{smear}{\left(b \right)}=T_{int,used} \max\left(1.0, \min\left(\frac{T_{snap}}{T_{int,used}}, \left\lfloor{\frac{\epsilon_{f approx} \lambda}{\Omega_{E} T_{int,used} \Theta_{fov} b}}\right\rfloor\right)\right)\end{aligned}$
Number of Channels at the Frequency Smearing Limit, Nf,smear¶
In [23]:
show_sympy(Symbol("Theta_PSF"), tp3.Theta_beam)
show_sympy(Symbol("Q_bw"), tp3.Qbw)
show_sympy(Function("N_f,smear")(b), tp2.Nf_no_smear_backward.term)
$\displaystyle \begin{aligned}\Theta_{PSF}=\frac{1.5 \lambda}{B_{max} \sqrt{Q_{subband}}}\end{aligned}$
$\displaystyle \begin{aligned}Q_{bw}=\frac{0.600124986981879}{\sqrt{1 - \frac{1}{amp_{f max}}}}\end{aligned}$
$\displaystyle \begin{aligned}N_{f,smear}{\left(b \right)}=\frac{\log{\left(\frac{\lambda_{max}}{\lambda_{min}} \right)}}{\log{\left(1 + \frac{1.5 \lambda}{Q_{bw} \Theta_{fov} b} \right)}}\end{aligned}$
Image Size¶
In [24]:
from numpy import pi
# Replace wl_sb_max and pi by symbol in formula
show_sympy(Symbol("Theta_FoV"), tp3.Theta_fov_total /tp3.wl_sb_max*tp.wl_sb_max * pi/Symbol("pi"))
$\displaystyle \begin{aligned}\Theta_{FoV}=\frac{7.66 Q_{fov} \lambda_{sb max}}{D_{s} \pi}\end{aligned}$
Image Plane Pixel Size¶
In [25]:
# Replace wl_sb_min by symbol in formula
show_sympy(Symbol("Theta_pix"), tp3.Theta_pix /tp3.wl_sb_min*tp.wl_sb_min)
$\displaystyle \begin{aligned}\Theta_{pix}=\frac{0.75 \lambda_{sb min}}{B_{max} Q_{pix}}\end{aligned}$
Number of Pixels on Image or Grid Size¶
In [26]:
# All replacements from above
show_sympy(Symbol("N_pix"), tp3.Npix_linear_fov_total
/tp3.wl_sb_max*tp.wl_sb_max *tp3.wl_sb_min/tp.wl_sb_min *pi/Symbol("pi"))
$\displaystyle \begin{aligned}N_{pix}=\frac{10.2133333333333 B_{max} Q_{fov} Q_{pix} \lambda_{sb max}}{D_{s} \lambda_{sb min} \pi}\end{aligned}$
Update Frequency Scale for Convolution Kernels¶
In [27]:
show_sympy(Function("N_f,kernel")(b), tp2.Nf_gcf_backward.term)
$\displaystyle \begin{aligned}N_{f,kernel}{\left(b \right)}=\max\left(N_{f,min}, \min\left(N_{f,max}, \frac{\log{\left(\frac{\lambda_{max}}{\lambda_{min}} \right)}}{\log{\left(1.0 + \frac{\epsilon_{f approx}}{Q_{kernel} \left(N_{aa}^{2} + 1.0 \Theta_{fov}^{2} \left(\Theta_{fov}^{2} \min\left(\Delta_{W stack}, \Delta_{W max}{\left(b \right)}\right)^{2} + \frac{0.636619772367581 \Theta_{fov} \min\left(\Delta_{W stack}, \Delta_{W max}{\left(b \right)}\right)^{1.5}}{\epsilon_{w}}\right)\right)^{0.5}} \right)}}\right)\right)\end{aligned}$
Size of the w-kernel¶
In [28]:
show_sympy(1/(2*Symbol("pi")),1/2/pi)
show_sympy(Function("N_GW")(b), tp2.Ngw_backward.term)
show_sympy(tp.DeltaW_max.term, tp3.DeltaW_max.term)
$\displaystyle \begin{aligned}\frac{1}{2 \pi}=0.159154943091895\end{aligned}$
$\displaystyle \begin{aligned}N_{GW}{\left(b \right)}=1.0 \Theta_{fov} \sqrt{\Theta_{fov}^{2} \min\left(\Delta_{W stack}, \Delta_{W max}{\left(b \right)}\right)^{2} + \frac{0.636619772367581 \Theta_{fov} \min\left(\Delta_{W stack}, \Delta_{W max}{\left(b \right)}\right)^{1.5}}{\epsilon_{w}}}\end{aligned}$
$\displaystyle \begin{aligned}\Delta_{W max}{\left(b \right)}=Q_{w} \max\left(\frac{0.125 b^{2}}{R_{Earth} \lambda}, \frac{0.5 b \sin{\left(\Omega_{E} T_{snap} \right)}}{\lambda}\right)\end{aligned}$
Imaging Pipeline Geometry Assumptions¶
In [29]:
show_sympy(Function("\Delta w_min")(b), tp3.DeltaW_Earth.term)
<>:1: SyntaxWarning: invalid escape sequence '\D' <>:1: SyntaxWarning: invalid escape sequence '\D' /tmp/ipykernel_175/535096630.py:1: SyntaxWarning: invalid escape sequence '\D' show_sympy(Function("\Delta w_min")(b), tp3.DeltaW_Earth.term)
$\displaystyle \begin{aligned}\Delta w_{min}{\left(b \right)}=\frac{0.125 b^{2}}{R_{Earth} \lambda}\end{aligned}$