SITCOMTN-139
WET-006: pairing vs stacking for wavefront estimation#
In this document we compare two approaches for wavefront estimation for the Rubin Observatory Active Optics System (AOS). The first (default) approach, dubbed “pairing”, is using individual donut pairs and averaging their final Zk estimate. The second approach (“stacking”) involves aggregating individual donut stamps prior to wavefront fitting and using their combined information (a stacked stamp) to obtain a single set of Zernike polynomials per stacked donut image. For an overview of pairing vs stacking, also see https://sitcomtn-085.lsst.io/.
Last verified to run: 09/19/2024
Versions:
lsst_distrib w_2024_32 (ext, cvmfs)
ts_wep v11.4.0
Setup on USDF#
We import all required packages:
import matplotlib.pyplot as plt
from astropy.visualization import ZScaleInterval
from lsst.daf import butler as dafButler
import matplotlib.pyplot as plt
import numpy as np
from lsst.ts.wep.utils import (
WfAlgorithmName,
getTaskInstrument,
getCameraFromButlerName
)
from lsst.ts.wep.image import Image
import stacking
Data processing#
The data was processed by running first ISR (lsstPipelineISRdoVar.yaml ):
description: run ISR on raws
instrument: lsst.obs.lsst.LsstCam
tasks:
isr:
class: lsst.ip.isr.isrTask.IsrTask
config:
connections.outputExposure: ‘postISRCCD’
doBias: False
doVariance: True
doLinearize: False
doCrosstalk: False
doDefect: False
doNanMasking: False
doInterpolate: False
doBrighterFatter: False
doDark: False
doFlat: False
doApplyGains: True
doFringe: False
doOverscan: True
python: OverscanCorrectionTask.ConfigClass.fitType = ‘MEDIAN’
Then WEP donut direct detection and cutout (lsstPipelineDirectCutoutOnly.yaml):
description: detect and cutout donuts
instrument: lsst.obs.lsst.LsstCam
tasks:
generateDonutDirectDetectTask:
class: lsst.ts.wep.task.generateDonutDirectDetectTask.GenerateDonutDirectDetectTask
config:
donutSelector.useCustomMagLimit: True
cutOutDonutsCwfsTask:
class: lsst.ts.wep.task.cutOutDonutsCwfsTask.CutOutDonutsCwfsTask
Then separately using the donut postage stamps to obtain Danish-based (lsstPipelineCalcDanishOnly.yaml):
description: calculate Zernikes with Danish
instrument: lsst.obs.lsst.LsstCam
tasks:
calcZernikesTask:
class: lsst.ts.wep.task.calcZernikesTask.CalcZernikesTask
config:
python: |
from lsst.ts.wep.task import EstimateZernikesDanishTask
config.estimateZernikes.retarget(EstimateZernikesDanishTask)
or TIE-based wavefront estimation (lsstPipelineCalcTieOnly.yaml):
description: calculate Zernikes with TIE
instrument: lsst.obs.lsst.LsstCam
tasks:
calcZernikesTask:
class: lsst.ts.wep.task.calcZernikesTask.CalcZernikesTask
config:
python: |
from lsst.ts.wep.task import EstimateZernikesTieTask
config.estimateZernikes.retarget(EstimateZernikesTieTask)
The succession of these was run with in a new repo where we first register the instrument, write calibrations, and ingest the raws:
import os
pathCwd = '/sdf/group/rubin/shared/scichris/DM-42718_WET-006'
butlerRootPath = os.path.join(pathCwd, 'newRepo')
outputImgDir = '/sdf/data/rubin/user/gmegias/projects/commissioning_sims/WET-001_lsstcam/state_*'
cmd1 = f"butler create {butlerRootPath}"
butlerInstName = 'Cam'
cmd2 = f"butler register-instrument {butlerRootPath} lsst.obs.lsst.Lsst{butlerInstName}"
cmd3 = f"butler write-curated-calibrations {butlerRootPath} lsst.obs.lsst.Lsst{butlerInstName}"
cmd4 = f"butler ingest-raws {butlerRootPath} {outputImgDir}/amp*"
print(cmd1, '\n')
print(cmd2, '\n')
print(cmd3, '\n')
print(cmd4, '\n')
butler create /sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo
butler register-instrument /sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo lsst.obs.lsst.LsstCam
butler write-curated-calibrations /sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo lsst.obs.lsst.LsstCam
butler ingest-raws /sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo /sdf/data/rubin/user/gmegias/projects/commissioning_sims/WET-001_lsstcam/state_*/amp*
Then we run ISR, and each WEP pipetasks separately so that they are in collections that can be reused. Each was run with bps on USDF, with the following prerequisites:
source /sdf/group/rubin/sw/w_latest/loadLSST.bash && setup lsst_distrib
source setup_aos.sh, which contains:
#!/usr/bin/env bash
export PATH_TO_TS_WEP=/sdf/data/rubin/user/scichris/WORK/aos_packages/ts_wep/
export PATH_TO_TS_OFC=/sdf/data/rubin/user/scichris/WORK/aos_packages/ts_ofc/
export PATH_TO_TS_IMSIM=/sdf/data/rubin/user/scichris/WORK/aos_packages/ts_imsim/
eups declare -r $PATH_TO_TS_WEP -t $USER --nolocks
eups declare -r $PATH_TO_TS_OFC -t $USER --nolocks
eups declare -r $PATH_TO_TS_IMSIM -t $USER --nolocks
setup ts_ofc -t $USER -t current
setup ts_wep -t $USER -t current
setup ts_imsim -t $USER -t current
export LSST_ALLOW_IMPLICIT_THREADS=True
Then
cd /sdf/group/rubin/shared/scichris/DM-42718_WET-006/
allocateNodes.py -v -n 10 -c 64 -m 60:00:00 -q roma,milano -g 1800 s3df --account rubin:commissioning
isrCollection = 'run_isr_all'
cutoutsCollection = 'all_states_direct_stamps'
tieCollection = 'all_states_direct_TIE'
danishCollection = 'all_states_direct_Danish'
isrYaml = 'lsstPipelineISRdoVar.yaml'
cutoutYaml = 'lsstPipelineDirectCutoutOnly.yaml'
danishYaml = 'lsstPipelineCalcDanishOnly.yaml'
tieYaml = 'lsstPipelineCalcTieOnly.yaml'
isrYamlPath = os.path.join(pathCwd, isrYaml)
cutoutYamlPath = os.path.join(pathCwd, cutoutYaml)
danishYamlPath = os.path.join(pathCwd, danishYaml)
tieYamlPath = os.path.join(pathCwd, tieYaml)
nCores = 50
cmdIsr = f'bps submit site_bps.yaml -b {butlerRootPath} -i LSSTCam/raw/all,LSSTCam/calib/unbounded \
-o {isrCollection} -p {isrYamlPath} -j {nCores}'
cmdCutout = f'bps submit site_bps.yaml -b {butlerRootPath} -i {isrCollection},LSSTCam/calib/unbounded \
-o {cutoutsCollection} -p {cutoutYamlPath} -j {nCores}'
cmdDanish= f'bps submit site_bps.yaml -b {butlerRootPath} -i {cutoutsCollection} \
-o {danishCollection} -p {danishYamlPath} -j {nCores}'
cmdTie= f'bps submit site_bps.yaml -b {butlerRootPath} -i {cutoutsCollection} \
-o {tieCollection} -p {tieYamlPath} -j {nCores}'
print(cmdIsr, '\n')
print(cmdCutout, '\n')
print(cmdDanish, '\n')
print(cmdTie, '\n')
bps submit site_bps.yaml -b /sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo -i LSSTCam/raw/all,LSSTCam/calib/unbounded -o run_isr_all -p /sdf/group/rubin/shared/scichris/DM-42718_WET-006/lsstPipelineISRdoVar.yaml -j 50
bps submit site_bps.yaml -b /sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo -i run_isr_all,LSSTCam/calib/unbounded -o all_states_direct_stamps -p /sdf/group/rubin/shared/scichris/DM-42718_WET-006/lsstPipelineDirectCutoutOnly.yaml -j 50
bps submit site_bps.yaml -b /sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo -i all_states_direct_stamps -o all_states_direct_Danish -p /sdf/group/rubin/shared/scichris/DM-42718_WET-006/lsstPipelineCalcDanishOnly.yaml -j 50
bps submit site_bps.yaml -b /sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo -i all_states_direct_stamps -o all_states_direct_TIE -p /sdf/group/rubin/shared/scichris/DM-42718_WET-006/lsstPipelineCalcTieOnly.yaml -j 50
Illustration of WET-001 states#
Note that each state contains randomly drawn modes:
butler = dafButler.Butler(butlerRootPath)
registry = butler.registry
datasetRefs = registry.queryDatasets('postISRCCD',collections=['run_isr_all'],
where=f"instrument='LSSTCam' ").expanded()
print(len(list(datasetRefs)))
dataIds = []
for ref in datasetRefs:
dataIds.append(ref.dataId)
800
# read in the images:
exposures = {}
for s in range(4):
exposures[s] = butler.get('postISRCCD',
dataId=dataIds[s],
collections=['run_isr_all']
)
# plot the example exposures
fig,ax = plt.subplots(1,4, figsize=(10,5))
zscale = ZScaleInterval()
i=0
for state, exposure in exposures.items():
d = exposure.image.array
#focusz = exposure.visitInfo.focusZ
vmin,vmax = zscale.get_limits(d)
mappable = ax[i].imshow(d.T, vmin=vmin, vmax=vmax, origin='lower')
ax[i].set_title(f'{exposure.visitInfo.instrumentLabel} {exposure.detector.getName()}, \nstate {state} ', )
ax[i].set_xlabel('x position [px]')
i += 1
ax[0].set_ylabel('y position [px]')
for i in range(1,4):
ax[i].set_yticklabels([])
fig.subplots_adjust(wspace=0.1)
Perform stacking#
Stacking can be performed with the code in the attached stacking.py file.
We find that for Danish there is a minimum background threshold, which means that it cannot be run entirely without backgroud
butlerRootPath = '/sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo'
output_collection = 'all_states_direct_stamps'
butler = dafButler.Butler(butlerRootPath)
dataId = {'instrument':'LSSTCam',
'detector':203,
'visit': 5024072302002
}
stampsIntra = butler.get('donutStampsIntra',
dataId=dataId,
collections=[output_collection]
)
fig,axs = plt.subplots(1,3,figsize=(12,4))
ax = np.ravel(axs)
i=0
use_mask = True
add_bkgnd = True
# set the file name title
string = 'no'
if use_mask: string = 'use'
bkgnd = 'no'
if add_bkgnd: bkgnd='with'
stackedIntra = stack_donut_wep_im_refactor(stampsIntra, N =len(stampsIntra),
pixel_stack = 'nanmean', use_mask = use_mask,
after_avg_fill_with_bkgnd=add_bkgnd
)
ax[i].imshow(stackedIntra['wep_im'].image, origin='lower',)
ax[i].set_title(f'{string} mask, add mean background')
i = 1
use_mask = True
add_bkgnd = False
# set the file name title
string = 'no'
if use_mask: string = 'use'
bkgnd = 'no'
if add_bkgnd: bkgnd='with'
stackedIntra = stack_donut_wep_im_refactor(stampsIntra, N =len(stampsIntra),
pixel_stack = 'nanmean', use_mask = use_mask,
after_avg_fill_with_bkgnd=add_bkgnd
)
ax[i].imshow(stackedIntra['wep_im'].image, origin='lower',)
ax[i].set_title(f'{string} mask, {bkgnd} background')
i = 2
use_mask = False
#add_bkgnd = False
# set the file name title
string = 'no'
if use_mask: string = 'use'
bkgnd = 'no'
if add_bkgnd: bkgnd='with'
stackedIntra = stack_donut_wep_im_refactor(stampsIntra, N =len(stampsIntra),
pixel_stack = 'nanmean', use_mask = use_mask,
after_avg_fill_with_bkgnd=add_bkgnd
)
ax[i].imshow(stackedIntra['wep_im'].image, origin='lower',)
ax[i].set_title(f'{string} mask: mean of image stack')
for i in range(3):
ax[i].set_xticks([])
ax[i].set_yticks([])
for i in [0,2]:
ax[i].text(50,200,'TIE $\checkmark $ Danish $\checkmark $', color='green', fontsize=13)
ax[1].text(50,200, 'TIE $\checkmark $' , color='green', fontsize=13)
ax[1].text(85,200, 'Danish X' , color='red', fontsize=13)
fig.suptitle('state 2, detector R44', fontsize=15)
/tmp/ipykernel_44053/2159255477.py:81: RuntimeWarning: Mean of empty slice
donutStackedArray = np.nanmean(image_stack, axis=0)
Text(0.5, 0.98, 'state 2, detector R44')
Illustrate the results#
Read all pre-computed stacking results into dictionary
summary = {}
use_mask = True
method = 'tie'
stacking_results = {}
# read the results for TIE, use mask, no bkgnd
fname = 'wep_direct_stacking_use_mask_tie_no_bkgnd.npy'
stacking_results[0] = {'use_mask':True, 'mask_string':'use',
'method':'tie',
'method_string': f'stacking_tie_use_mask',
'add_bkgnd':False, 'bkgnd_string': 'no',
'fname' :fname,
'results':np.load(fname, allow_pickle=True).item()
}
# read the results for {TIE, Danish} x {use / no mask} , with background added (if using mask)
i=1
for use_mask in [True, False]:
for method in ['tie', 'danish']:
add_bkgnd = True
# set the file name title
mask_string = 'no'
if use_mask: mask_string = 'use'
bkgnd_string = 'no'
if add_bkgnd: bkgnd_string = 'with'
fname = f'wep_direct_stacking_{mask_string}_mask_{method}_{bkgnd_string}_bkgnd.npy'
results = np.load(fname, allow_pickle=True).item()
method_string = f'stacking_{method}_{mask_string}_mask'
if mask_string == 'use' :
method_string += f'_{bkgnd_string}_bkgnd'
stacking_results[i] = {'use_mask':use_mask,
'mask_string':mask_string,
'method':method,
'method_string':method_string,
'add_bkgnd':add_bkgnd,
'bkgnd_string':bkgnd_string,
'fname':fname,
'results':results
}
i += 1
Read the pairing results from butler:
butlerRootPath = '/sdf/group/rubin/shared/scichris/DM-42718_WET-006/newRepo'
butler = dafButler.Butler(butlerRootPath)
registry = butler.registry
output_collection = 'all_states_direct_TIE'
datasetRefs = registry.queryDatasets('donutStampsExtra', collections=[output_collection],
where=f"instrument='LSSTCam' and detector.purpose='WAVEFRONT' ").expanded()
# where=f"instrument='LSST{butlerInstName}' and detector in (199) ").expanded()
print(len(list(datasetRefs)))
refs = []
for ref in datasetRefs:
refs.append(ref)
print(len(refs))
#stackedIntra
pairing_results = {'tie':{}, 'danish':{}}
for method in pairing_results.keys():
pairing_results[method] = {}
for state in range(1,101):
pairing_results[method][state] = {}
for ref in refs:
#print(ref.dataId)
# read in the results of each method
for method in pairing_results.keys():
if method == 'tie':
coll_method = 'TIE'
else:
coll_method = 'Danish'
coll = f'all_states_direct_{coll_method}'
state = int(str(ref.dataId.visit.id)[-3:])
raft = ref.dataId.detector.raft
pairing_results[method][state][raft] = butler.get('zernikeEstimateRaw',
dataId=ref.dataId,
collections=[coll])
400
400
Also load OPDs that are used as a reference:
# load all OPDs to a dict
from astropy.io import fits
all_opd = {}
for n in range(1,101):#101):
opdDir = f'/sdf/data/rubin/user/gmegias/projects/commissioning_sims/WET-001_lsstcam/state_{n}'
hdul = fits.open(os.path.join(opdDir,'opd.fits'))
opds = {}
for i in range(len(hdul)):
opd_zks_1_28 = []
for key,value in hdul[i].header.items():
if key.startswith('AZ'):
#print(key.split('_')[1])
opd_zks_1_28.append(value)
opds[i] = opd_zks_1_28
all_opd[n] = opds
opd_raft_to_id = {'R00':0, 'R04':1, 'R40':2, 'R44':3}
Compare TIE to Danish for pairing only#
import matplotlib.pyplot as plt
import numpy as np
rafts = ['R04', 'R44', 'R00','R40']
state=5
fig,axs = plt.subplots(2,2,figsize=(16,8))
ax = np.ravel(axs)
i=0
# this plots just the pairing results
colors = {'danish':'green', 'tie':'orange'}
for raft in rafts:
ax[i].set_title(raft )
opdId = opd_raft_to_id[raft]
opd_zk = 0.001*np.array(all_opd[state][opdId][3:22])
for method in pairing_results.keys():
# plot pairing raw zks for danish and tie
zkRaw = pairing_results[method][state][raft]
avg_zk = np.mean(zkRaw, axis=0)
rms_diff_avg = np.sqrt(np.mean(np.square(avg_zk-opd_zk)))
ax[i].plot(np.arange(4,23), avg_zk-opd_zk, marker='d', label=f'{method}')
ax[i].set_ylabel('Delta Zk mode [microns]')
ax[i].set_xlabel('Zk mode')
ax[i].set_xticks([4,6,8,10,12,14,16,18,20,22])
ax[i].axhline(0,ls='--', c='red')
i+=1
ax[0].legend()#bbox_to_anchor=[1.2,0.8])
fig.subplots_adjust(hspace=0.3)
fig.suptitle(f'WET-006 lsstCam, state {state}, \n pairing only (average of zernikeEstimateRaw) : TIE vs Danish')
Text(0.5, 0.98, 'WET-006 lsstCam, state 5, \n pairing only (average of zernikeEstimateRaw) : TIE vs Danish')
Add to that the stacking results for that particular mode, and calculate the RMS difference:
import matplotlib.pyplot as plt
import numpy as np
rafts = ['R04', 'R44', 'R00','R40']
state=5
fig,axs = plt.subplots(2,2,figsize=(16,8))
ax = np.ravel(axs)
i=0
# this plots the two sets of results
colors = {'danish':'green', 'tie':'orange'}
for raft in rafts:
ax[i].set_title(raft )
# get OPD for that raft
opdId = opd_raft_to_id[raft]
opd_zk = 0.001*np.array(all_opd[state][opdId][3:22])
for method in pairing_results.keys():
zkRaw = pairing_results[method][state][raft]
avg_zk = np.mean(zkRaw, axis=0)
rms_diff_avg = np.sqrt(np.mean(np.square(avg_zk-opd_zk)))
ax[i].set_title(f'{raft} '+r'$\Delta$ rms '+\
f'{method} (stack: {rms_diff_stack:.3f}, ' +\
f'avg: {rms_diff_avg:.3f} )')
ax[i].plot(np.arange(4,23), avg_zk-opd_zk, marker='d',
label=f'{method} (avg - OPD)')
for method in stacking_results.keys():
stacked_zk = stacking_results[method]['results'][state][raft]
label = stacking_results[method]['method_string']
rms_diff_stack = np.sqrt(np.mean(np.square(stacked_zk-opd_zk)))
ax[i].plot(np.arange(4,23), stacked_zk-opd_zk, marker='d',
label=f'{label} (stacked - OPD)')
ax[i].set_ylabel('Delta Zk mode [microns]')
ax[i].set_xlabel('Zk mode')
ax[i].set_xticks([4,6,8,10,12,14,16,18,20,22])
ax[i].axhline(0,ls='--', c='red')
i+=1
ax[1].legend()
fig.subplots_adjust(hspace=0.3)
fig.suptitle(f'WET-006 lsstCam, state {state}, RMS difference (fit - OPD) ')
Text(0.5, 0.98, 'WET-006 lsstCam, state 5, RMS difference (fit - OPD) ')
Marginalize this information across all states:
rafts = ['R04', 'R44', 'R00','R40']
rmss={}
# add stacking results
for j in stacking_results.keys():
results = stacking_results[j]['results']
mask_string = stacking_results[j]['mask_string']
bkgnd_string = stacking_results[j]['bkgnd_string']
method = stacking_results[j]['method']
method_string = f'stacking_{method}_{mask_string}_mask'
if mask_string == 'use' :
method_string += f'_{bkgnd_string}_bkgnd'
print(method_string)
# else:
# i=0
#method+='_stacking_'
rmss[method_string] = {}
for raft in rafts:
rmss[method_string][raft] = []
opdId = opd_raft_to_id[raft]
for state in range(1,101):
opd_zk = 0.001*np.array(all_opd[state][opdId][3:22])
stacked_zk = results[state][raft]
rms_diff = np.sqrt(np.mean(np.square(stacked_zk-opd_zk)))
rmss[method_string][raft].append(rms_diff)
# add pairing results
for algo in pairing_results.keys():
method_string = f'pairing_{algo}'
rmss[method_string] = {}
print(method_string)
for raft in rafts:
rmss[method_string][raft] = []
opdId = opd_raft_to_id[raft]
for state in range(1,101):
opd_zk = 0.001*np.array(all_opd[state][opdId][3:22])
avg_zk = np.mean(pairing_results[algo][state][raft], axis=0)
rms_diff = np.sqrt(np.mean(np.square(avg_zk-opd_zk)))
rmss[method_string][raft].append(rms_diff)
stacking_tie_use_mask_no_bkgnd
stacking_tie_use_mask_with_bkgnd
stacking_danish_use_mask_with_bkgnd
stacking_tie_no_mask
stacking_danish_no_mask
pairing_tie
pairing_danish
rafts = ['R04', 'R44', 'R00','R40']
fig,axs = plt.subplots(2,2,figsize=(16,10))
ax = np.ravel(axs)
for method in rmss.keys():
i=0
for raft in rafts:
ax[i].plot(np.arange(1,101), rmss[method][raft], marker='s', label=f'{method}')
i+=1
fig.text(0.5,0.06,'state #')
fig.text(0.06,0.5,r'RMS ($\Delta$(zk-opd) ) ', rotation=90)
fig.suptitle('pairing vs all stacking methods', fontsize=16)
ax[3].legend(bbox_to_anchor=[1.05,1.2])
<matplotlib.legend.Legend at 0x7f97689ff0d0>
This doesn’t seem to prefer either method, and to make it easier to see the difference we calculate the mean RMS difference across all states per detector per method:
from matplotlib import colormaps as cmaps
from astropy.table import Table
means_arr = []
rafts_arr = []
methods_arr = []
rafts = ['R04', 'R44', 'R00','R40']
for raft in rafts:
#means[raft] =
for method in rmss.keys():
mean_rms = np.median( rmss[method][raft])
#print(raft, mean_rms, f'{method}')
means_arr.append(mean_rms)
rafts_arr.append(raft)
methods_arr.append(method)
#print('\n')
mean_rms_table = Table(data=[means_arr, rafts_arr, methods_arr], names=['mean_rms', 'raft', 'method'])
for raft in rafts:
per_raft = mean_rms_table[mean_rms_table['raft'] == raft]# .sort('mean_rms')
per_raft.sort('mean_rms')
print(per_raft)
mean_rms raft method
------------------- ---- -----------------------------------
0.12245982433194912 R04 stacking_tie_use_mask_with_bkgnd
0.12287115478025104 R04 stacking_tie_no_mask
0.1272135180411988 R04 pairing_tie
0.1337388102052313 R04 stacking_danish_no_mask
0.1342364041823671 R04 stacking_danish_use_mask_with_bkgnd
0.13711868256738124 R04 pairing_danish
0.16788752144434005 R04 stacking_tie_use_mask_no_bkgnd
mean_rms raft method
------------------- ---- -----------------------------------
0.13367654248935212 R44 stacking_danish_use_mask_with_bkgnd
0.13394723071818743 R44 stacking_tie_use_mask_with_bkgnd
0.15050260405928126 R44 pairing_danish
0.1521178698309197 R44 stacking_danish_no_mask
0.1597151828776096 R44 stacking_tie_use_mask_no_bkgnd
0.1836104280153719 R44 pairing_tie
0.18693449833802006 R44 stacking_tie_no_mask
mean_rms raft method
------------------- ---- -----------------------------------
0.09211418516055045 R00 pairing_danish
0.1066278700894988 R00 stacking_danish_use_mask_with_bkgnd
0.11628452906679443 R00 stacking_danish_no_mask
0.11727095036248558 R00 stacking_tie_use_mask_no_bkgnd
0.13644630985049472 R00 stacking_tie_use_mask_with_bkgnd
0.15912864847592118 R00 pairing_tie
0.2085669551915312 R00 stacking_tie_no_mask
mean_rms raft method
------------------- ---- -----------------------------------
0.09740636512101784 R40 stacking_danish_no_mask
0.10266305805150648 R40 stacking_danish_use_mask_with_bkgnd
0.10726847179842722 R40 pairing_danish
0.11828681239339955 R40 pairing_tie
0.14696597370236864 R40 stacking_tie_use_mask_with_bkgnd
0.20589407061432713 R40 stacking_tie_no_mask
0.2073340524989511 R40 stacking_tie_use_mask_no_bkgnd
We can also marginalize across detectors, summarizing the RMS difference between each method taking the average per four sensors:
mean_per_state = {}
for method in rmss.keys():
all_raft_data = []
for raft in rafts:
all_raft_data.append( rmss[method][raft])
mean_per_state[method] = np.mean(all_raft_data,axis=0)
Plot per state:
import matplotlib.pyplot as plt
from matplotlib import colormaps as cmaps
rafts = ['R04', 'R44', 'R00','R40']
cmap = cmaps['tab10']
fig,ax = plt.subplots(1,1,figsize=(8,4))
j=0
for method in mean_per_state.keys():
ax.plot(np.arange(1,101), mean_per_state[method], marker='o', label=f'{method}',
c=cmap(j))
j+=1
ax.set_xlabel('state #')
ax.set_ylabel(r'mean RMS ($\Delta$(zk-opd) across rafts) ',)
ax.set_title('pairing vs all stacking methods', fontsize=16)
ax.legend(bbox_to_anchor=[1.55,0.75])
<matplotlib.legend.Legend at 0x7f97656f3650>
Show just the two pairing methods:
import matplotlib.pyplot as plt
rafts = ['R04', 'R44', 'R00','R40']
fig,ax = plt.subplots(1,1,figsize=(8,4))
cmap = cmaps['tab10']
j=0
for method in mean_per_state.keys():
if method.startswith('pairing'):
ax.plot(np.arange(1,101), mean_per_state[method], marker='o', label=f'{method}', c=cmap(j))
j +=1
ax.set_xlabel('state #')
ax.set_ylabel(r'mean RMS ($\Delta$(zk-opd) across rafts ',)
ax.set_title('pairing methods', fontsize=16)
ax.legend()#bbox_to_anchor=[1.55,0.75])
<matplotlib.legend.Legend at 0x7f97685b5c50>
Danish appears to overall produce results closer to the OPD.
Most stacking results are not as good as pairing for Danish:
import matplotlib.pyplot as plt
rafts = ['R04', 'R44', 'R00','R40']
fig,ax = plt.subplots(1,1,figsize=(8,4))
cmap = cmaps['tab10']
j=0
for method in mean_per_state.keys():
if 'danish' in method:
ax.plot(np.arange(1,101), mean_per_state[method], marker='o', label=f'{method}', c=cmap(j))
j +=1
ax.set_xlabel('state #')
ax.set_ylabel(r'mean RMS ($\Delta$(zk-opd) across rafts ',)
ax.set_title('pairing vs stacking methods for Danish', fontsize=16)
ax.legend()#bbox_to_anchor=[1.55,0.75])
<matplotlib.legend.Legend at 0x7f97684060d0>
And for TIE stacking also does provide statistically significant improvement over pairing, the RMS floor being around 0.09 for both approaches
import matplotlib.pyplot as plt
rafts = ['R04', 'R44', 'R00','R40']
fig,ax = plt.subplots(1,1,figsize=(8,4))
cmap = cmaps['tab10']
j=0
for method in mean_per_state.keys():
if 'tie' in method:
ax.plot(np.arange(1,101), mean_per_state[method], marker='o', label=f'{method}', c=cmap(j))
j +=1
ax.set_xlabel('state #')
ax.set_ylabel(r'mean RMS ($\Delta$(zk-opd) across rafts ',)
ax.set_title('pairing vs stacking methods for TIE', fontsize=16)
ax.legend()#bbox_to_anchor=[1.55,0.75])|
<matplotlib.legend.Legend at 0x7f9768404150>
In summary, based on the WET-001 simulation we do not see that there would be a clear advantage in introducing stacking over using the existing approach of pairing.