Skip to content

Commit 3369a53

Browse files
ajeldoradoajeldorado
authored
Merge pull request #106 from ajeldorado/algodiff_aj
Pairwise probing updates
2 parents 43b592a + 6fc6452 commit 3369a53

3 files changed

Lines changed: 63 additions & 36 deletions

File tree

‎falco/est.py‎

Lines changed: 48 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -240,7 +240,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
240240
Itr = ev.Itr
241241
whichDM = mp.est.probe.whichDM
242242

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

294294
# Store the initial DM commands
295295
if np.any(mp.dm_ind == 1):
296-
DM1Vnom = mp.dm1.V
296+
DM1Vnom = mp.dm1.V.copy()
297297
else:
298298
DM1Vnom = np.zeros((mp.dm1.Nact, mp.dm1.Nact))
299299

300300
if np.any(mp.dm_ind == 2):
301-
DM2Vnom = mp.dm2.V
301+
DM2Vnom = mp.dm2.V.copy()
302302
else:
303303
DM2Vnom = np.zeros((mp.dm2.Nact, mp.dm2.Nact))
304304

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

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

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

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

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

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

545547
# For unprobed field based on model:
546548
if np.any(mp.dm_ind == 1):
547-
mp.dm1.V = DM1Vnom
549+
mp.dm1.V = DM1Vnom.copy()
548550
if np.any(mp.dm_ind == 2):
549-
mp.dm2.V = DM2Vnom
551+
mp.dm2.V = DM2Vnom.copy()
550552

551553
E0 = falco.model.compact(mp, modvar)
552554
E0vec = E0[mp.Fend.corr.maskBool]
@@ -578,8 +580,8 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
578580
dphdm = np.zeros((mp.Fend.corr.Npix, Npairs)) # phases
579581
for iProbe in range(Npairs):
580582
dphdm[:, iProbe] = np.arctan2(
581-
np.imag(dEplus[:, iProbe]) - np.imag(dEminus[:, iProbe]),
582-
np.real(dEplus[:, iProbe]) - np.real(dEminus[:, iProbe]))
583+
np.imag(dEplus[:, iProbe] - dEminus[:, iProbe]),
584+
np.real(dEplus[:, iProbe] - dEminus[:, iProbe]))
583585

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

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

597600
if mp.est.flagUseJac:
598601
dE = dEplus[ipix, :].T
599602
H = np.array([np.real(dE), np.imag(dE)])
603+
Epix = np.linalg.pinv(H) @ zAll[:, ipix] # Batch processing
600604
else:
605+
allProbeInds = np.arange(Npairs, dtype=int)
606+
goodProbeInds = allProbeInds[amp[ipix, :] > 0]
607+
NpairsGood = len(goodProbeInds)
608+
601609
H = np.zeros([Npairs, 2]) # Observation matrix
602-
# Leave Eest for a pixel as zero if any probe amp is 0
603-
if isnonzero[ipix] == 1:
604-
for iProbe in range(Npairs):
605-
H[iProbe, :] = amp[ipix, iProbe] * \
606-
np.array([np.cos(dphdm[ipix, iProbe]),
607-
np.sin(dphdm[ipix, iProbe])])
608-
else:
610+
611+
# If <2 probe pairs had good measurements, can't do pinv. Leave Eest as zero.
612+
if NpairsGood < 2:
613+
zerosCounter = zerosCounter + 1
614+
Epix = np.zeros((2, 1))
609615
zerosCounter += 1
610616

611-
Epix = np.linalg.pinv(H) @ zAll[:, ipix] # Batch processing
617+
# Otherwise, use the 2+ good probe pair measurements for that pixel
618+
else:
619+
for iProbe in range(Npairs):
620+
H[iProbe, :] = amp[ipix, iProbe] * np.array([
621+
np.cos(dphdm[ipix, iProbe]), np.sin(dphdm[ipix, iProbe])])
622+
623+
Epix = np.linalg.pinv(H[goodProbeInds, :]) @ zAll[goodProbeInds, ipix] # Batch processing
624+
625+
# Record how many pixels didn't use all the probe pairs
626+
if NpairsGood < Npairs:
627+
partialCounter = partialCounter + 1
628+
629+
# # Leave Eest for a pixel as zero if any probe amp is 0
630+
# if isnonzero[ipix] == 1:
631+
# for iProbe in range(Npairs):
632+
# H[iProbe, :] = amp[ipix, iProbe] * \
633+
# np.array([np.cos(dphdm[ipix, iProbe]),
634+
# np.sin(dphdm[ipix, iProbe])])
635+
612636
Eest[ipix] = Epix[0] + 1j*Epix[1]
613637

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

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

620647
# Initialize the state and state covariance estimates for the
621648
# Kalman filter. The state is the real and imag parts of the

‎falco/plot.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -317,7 +317,7 @@ def pairwise_probes(mp, ev, dDMVplus, ampSq2Dcube, iSubband):
317317
for row in range(Npairs):
318318
ax = axs[row, col]
319319
ax.set_title(titles[col])
320-
ax.invert_yaxis()
320+
# ax.invert_yaxis()
321321
ax.set_box_aspect(1)
322322
if row != Npairs-1:
323323
ax.tick_params(labelbottom=False, bottom=False)

‎tests/functional/test_wfsc_flc.py‎

Lines changed: 14 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -23,32 +23,32 @@ def test_wfsc_flc():
2323
out = falco.setup.flesh_out_workspace(mp)
2424
falco.wfsc.loop(mp, out)
2525

26-
# print(out.IrawCorrHist[-1])
27-
# print(out.IestScoreHist[-1])
28-
# print(out.IincoCorrHist[-1])
29-
# print(out.complexProjection[1, 0])
30-
# print(out.dm1.Spv[-1])
31-
# print(out.thput[-1])
32-
# print(out.log10regHist)
26+
print(out.IrawCorrHist[-1])
27+
print(out.IestScoreHist[-1])
28+
print(out.IincoCorrHist[-1])
29+
print(out.complexProjection[1, 0])
30+
print(out.dm1.Spv[-1])
31+
print(out.thput[-1])
32+
print(out.log10regHist)
3333

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

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

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

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

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

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

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

0 commit comments

Comments
 (0)