Skip to content

Pairwise probing updates #106

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ajeldorado merged 3 commits into from
Jul 18, 2025
Merged

Pairwise probing updates #106

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
14eaefa
Add estimator logic to keep a pixel's probes if at least 2 probe pair…
ajeldorado Jul 18, 2025
04aa446
Updated functional regression test now that pairwise probing throws o…
ajeldorado Jul 18, 2025
6fc6452
Set plotting to off in functional test
ajeldorado Jul 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 48 additions & 21 deletions falco/est.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -240,7 +240,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
Itr = ev.Itr
whichDM = mp.est.probe.whichDM

# # Reset "ev" if not using a Kalman filter.
# # Reset "ev" if not using a Kalman filter.
# if mp.estimator != 'pwp-kf':
# # ev = falco.config.Object()
# ev.dm1 = falco.config.Object()
Expand Down Expand Up @@ -293,12 +293,12 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):

# Store the initial DM commands
if np.any(mp.dm_ind == 1):
DM1Vnom = mp.dm1.V
DM1Vnom = mp.dm1.V.copy()
else:
DM1Vnom = np.zeros((mp.dm1.Nact, mp.dm1.Nact))

if np.any(mp.dm_ind == 2):
DM2Vnom = mp.dm2.V
DM2Vnom = mp.dm2.V.copy()
else:
DM2Vnom = np.zeros((mp.dm2.Nact, mp.dm2.Nact))

Expand All @@ -323,7 +323,9 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
for k in range(Npairs-1):
probePhaseVec = np.append(probePhaseVec, probePhaseVec[-1]-(Npairs-1))
probePhaseVec = np.append(probePhaseVec, probePhaseVec[-1]+Npairs)
# print(f'probePhaseVec = {probePhaseVec}')
probePhaseVec = probePhaseVec*np.pi/(Npairs)
print(f'probePhaseVec = {probePhaseVec}')

# Only for square probe region centered on the star
if mp.estimator in ('pairwise', 'pairwise-square', 'pwp-bp-square'):
Expand Down Expand Up @@ -378,8 +380,8 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):

# Measure current contrast level average
# Reset DM commands to the unprobed state:
mp.dm1.V = DM1Vnom
mp.dm2.V = DM2Vnom
mp.dm1.V = DM1Vnom.copy()
mp.dm2.V = DM2Vnom.copy()
# Separate out image values at DH pixels and delta DM voltage settings
Iplus = np.zeros((mp.Fend.corr.Npix, Npairs))
Iminus = np.zeros((mp.Fend.corr.Npix, Npairs))
Expand All @@ -391,7 +393,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
# Compute probe shapes and take probed images:

# Take initial, unprobed image (for unprobed DM settings).
whichImage = 1
whichImage = 0
I0 = falco.imaging.get_sbp_image(mp, iSubband)
I0vec = I0[mp.Fend.corr.maskBool] # Vectorize the dark hole

Expand All @@ -402,9 +404,9 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
# Store values for first image and its DM commands
ev.imageArray[:, :, whichImage, iSubband] = I0
if np.any(mp.dm_ind == 1):
ev.dm1.Vall[:, :, whichImage, iSubband] = mp.dm1.V
ev.dm1.Vall[:, :, whichImage, iSubband] = mp.dm1.V.copy()
if np.any(mp.dm_ind == 2):
ev.dm2.Vall[:, :, whichImage, iSubband] = mp.dm2.V
ev.dm2.Vall[:, :, whichImage, iSubband] = mp.dm2.V.copy()

# Compute the average Inorm in the scoring and correction regions
ev.corr = falco.config.Object()
Expand Down Expand Up @@ -502,7 +504,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
ampSq[ampSq < 0] = 0 # If probe amplitude is zero, set amp = 0
amp = np.sqrt(ampSq) # E-field amplitudes, dimensions: [mp.Fend.corr.Npix, Npairs]
isnonzero = np.all(amp, 1)
zAll = ((Iplus-Iminus)/4).T # Measurement vector, dimensions: [Npairs,mp.Fend.corr.Npix]
zAll = ((Iplus-Iminus)/4).T # Measurement vector, dimensions: [Npairs, mp.Fend.corr.Npix]
ampSq2Dcube = np.zeros((mp.Fend.Neta, mp.Fend.Nxi, mp.est.probe.Npairs))
for iProbe in range(Npairs): # Display the actual probe intensity
ampSq2D = np.zeros((mp.Fend.Neta, mp.Fend.Nxi))
Expand Down Expand Up @@ -544,9 +546,9 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):

# For unprobed field based on model:
if np.any(mp.dm_ind == 1):
mp.dm1.V = DM1Vnom
mp.dm1.V = DM1Vnom.copy()
if np.any(mp.dm_ind == 2):
mp.dm2.V = DM2Vnom
mp.dm2.V = DM2Vnom.copy()

E0 = falco.model.compact(mp, modvar)
E0vec = E0[mp.Fend.corr.maskBool]
Expand Down Expand Up @@ -578,8 +580,8 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
dphdm = np.zeros((mp.Fend.corr.Npix, Npairs)) # phases
for iProbe in range(Npairs):
dphdm[:, iProbe] = np.arctan2(
np.imag(dEplus[:, iProbe]) - np.imag(dEminus[:, iProbe]),
np.real(dEplus[:, iProbe]) - np.real(dEminus[:, iProbe]))
np.imag(dEplus[:, iProbe] - dEminus[:, iProbe]),
np.real(dEplus[:, iProbe] - dEminus[:, iProbe]))

# Batch process the measurements to estimate the electric field in the
# dark hole. Done pixel by pixel.
Expand All @@ -592,30 +594,55 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):

Eest = np.zeros((mp.Fend.corr.Npix,), dtype=complex)
zerosCounter = 0 # number of zeroed-out dark hole pixels
partialCounter = 0 # number of pixels with >=2 but <Npairs usable probe pairs.
for ipix in range(mp.Fend.corr.Npix):

if mp.est.flagUseJac:
dE = dEplus[ipix, :].T
H = np.array([np.real(dE), np.imag(dE)])
Epix = np.linalg.pinv(H) @ zAll[:, ipix] # Batch processing
else:
allProbeInds = np.arange(Npairs, dtype=int)
goodProbeInds = allProbeInds[amp[ipix, :] > 0]
NpairsGood = len(goodProbeInds)

H = np.zeros([Npairs, 2]) # Observation matrix
# Leave Eest for a pixel as zero if any probe amp is 0
if isnonzero[ipix] == 1:
for iProbe in range(Npairs):
H[iProbe, :] = amp[ipix, iProbe] * \
np.array([np.cos(dphdm[ipix, iProbe]),
np.sin(dphdm[ipix, iProbe])])
else:

# If <2 probe pairs had good measurements, can't do pinv. Leave Eest as zero.
if NpairsGood < 2:
zerosCounter = zerosCounter + 1
Epix = np.zeros((2, 1))
zerosCounter += 1

Epix = np.linalg.pinv(H) @ zAll[:, ipix] # Batch processing
# Otherwise, use the 2+ good probe pair measurements for that pixel
else:
for iProbe in range(Npairs):
H[iProbe, :] = amp[ipix, iProbe] * np.array([
np.cos(dphdm[ipix, iProbe]), np.sin(dphdm[ipix, iProbe])])

Epix = np.linalg.pinv(H[goodProbeInds, :]) @ zAll[goodProbeInds, ipix] # Batch processing

# Record how many pixels didn't use all the probe pairs
if NpairsGood < Npairs:
partialCounter = partialCounter + 1

# # Leave Eest for a pixel as zero if any probe amp is 0
# if isnonzero[ipix] == 1:
# for iProbe in range(Npairs):
# H[iProbe, :] = amp[ipix, iProbe] * \
# np.array([np.cos(dphdm[ipix, iProbe]),
# np.sin(dphdm[ipix, iProbe])])

Eest[ipix] = Epix[0] + 1j*Epix[1]

# If estimate is too bright, the estimate was probably bad.
Eest[np.abs(Eest)**2 > mp.est.Ithreshold] = 0.0

print('%d of %d pixels were given zero probe amplitude.' %
(zerosCounter, mp.Fend.corr.Npix))
if Npairs > 2: # Only possible for Npairs > 2
print('%d of %d pixels threw out at least one probe pair.' %
(partialCounter, mp.Fend.corr.Npix))

# Initialize the state and state covariance estimates for the
# Kalman filter. The state is the real and imag parts of the
Expand Down
2 changes: 1 addition & 1 deletion falco/plot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -317,7 +317,7 @@ def pairwise_probes(mp, ev, dDMVplus, ampSq2Dcube, iSubband):
for row in range(Npairs):
ax = axs[row, col]
ax.set_title(titles[col])
ax.invert_yaxis()
# ax.invert_yaxis()
ax.set_box_aspect(1)
if row != Npairs-1:
ax.tick_params(labelbottom=False, bottom=False)
Expand Down
28 changes: 14 additions & 14 deletions tests/functional/test_wfsc_flc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,32 +23,32 @@ def test_wfsc_flc():
out = falco.setup.flesh_out_workspace(mp)
falco.wfsc.loop(mp, out)

# print(out.IrawCorrHist[-1])
# print(out.IestScoreHist[-1])
# print(out.IincoCorrHist[-1])
# print(out.complexProjection[1, 0])
# print(out.dm1.Spv[-1])
# print(out.thput[-1])
# print(out.log10regHist)
print(out.IrawCorrHist[-1])
print(out.IestScoreHist[-1])
print(out.IincoCorrHist[-1])
print(out.complexProjection[1, 0])
print(out.dm1.Spv[-1])
print(out.thput[-1])
print(out.log10regHist)

# Tests:
Iend = out.IrawCorrHist[-1] # 1.1157e-5 in matlab, 1.073e-05 in python
assert isclose(Iend, 1.073e-5, abs_tol=1e-6)
assert isclose(Iend, 8.12e-6, abs_tol=1e-6)

Iest = out.IestScoreHist[-1] # 8.15e-06 in matlab, 8.123e-6 in python
assert isclose(Iest, 8.123e-6, abs_tol=3e-7)
assert isclose(Iest, 1.53e-5, abs_tol=3e-7)

Iinco = out.IincoCorrHist[-1] # 1.28e-5 in matlab, 1.165e-5 in python
assert isclose(Iinco, 1.165e-5, abs_tol=3e-7)
# Iinco = out.IincoCorrHist[-1] # 1.28e-5 in matlab, 1.165e-5 in python
# assert isclose(Iinco, 1.165e-5, abs_tol=3e-7)

complexProj = out.complexProjection[1, 0] # 0.74 in matlab, 0.82 in python
assert isclose(complexProj, 0.82, abs_tol=2e-2)
assert isclose(complexProj, 0.88, abs_tol=0.02)

dm1pv = out.dm1.Spv[-1] # 5.6956e-08 in matlab, 5.75e-8 in python
assert isclose(dm1pv, 5.75e-8, abs_tol=1e-9)
assert isclose(dm1pv, 6.09e-8, abs_tol=1e-9)

thput = out.thput[-1] # 0.1493 in matlab, 0.1492 in python
assert isclose(thput, 0.1492, abs_tol=1e-3)
assert isclose(thput, 0.1486, abs_tol=1e-3)

assert np.allclose(out.log10regHist, np.array([-2, -2, -2]), rtol=1e-2)

Expand Down