Skip to content

Commit 5100c74

Browse files
ajeldoradoajeldorado
committed
Add different relay value checks into the LC Jacobian unit test.
1 parent fb874a9 commit 5100c74

2 files changed

Lines changed: 123 additions & 65 deletions

File tree

‎tests/functional/test_jacobian_lc.py‎

Lines changed: 122 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -5,75 +5,133 @@
55

66
import falco
77

8-
import config_wfsc_lc as CONFIG
8+
import config_wfsc_lc_quick as CONFIG
99

1010

1111
def test_jacobian_lc():
1212
"""Lyot Jacobian test using FFTs between DMs."""
13-
mp = deepcopy(CONFIG.mp)
14-
mp.runLabel = 'test_lc'
15-
16-
mp.path = falco.config.Object()
17-
LOCAL_PATH = os.path.dirname(os.path.abspath(__file__))
18-
mp.path.falco = os.path.dirname(os.path.dirname(LOCAL_PATH))
19-
_ = falco.setup.flesh_out_workspace(mp)
20-
21-
# Fast Jacobian calculation
22-
mp.dm1.V = np.zeros((mp.dm1.Nact, mp.dm1.Nact))
23-
mp.dm2.V = np.zeros((mp.dm2.Nact, mp.dm2.Nact))
24-
jacStruct = falco.model.jacobian(mp) # Get structure containing Jacobians
25-
26-
G1fastAll = jacStruct.G1
27-
G2fastAll = jacStruct.G2
28-
Nind = 20
29-
subinds = np.ix_(np.arange(3, 20*4, 4).astype(int))
30-
absG1sum = np.sum(np.abs(G1fastAll), axis=1)
31-
indG1 = np.nonzero(absG1sum > 1e-2*np.max(absG1sum))[0]
32-
indG1subset = indG1[subinds] # Take a 20-actuator subset
33-
absG2sum = np.sum(np.abs(G2fastAll), axis=1)
34-
indG2 = np.nonzero(absG2sum > 1e-2*np.max(absG2sum))[0]
35-
indG2subset = indG2[subinds] # Take a 20-actuator subset
36-
G1fast = np.squeeze(G1fastAll[:, indG1subset])
37-
G2fast = np.squeeze(G2fastAll[:, indG2subset])
38-
39-
# Compute Jacobian via differencing (slower)
40-
modvar = falco.config.ModelVariables()
41-
modvar.whichSource = 'star'
42-
modvar.sbpIndex = 0
43-
modvar.starIndex = 0
44-
Eunpoked2D = falco.model.compact(mp, modvar)
45-
Eunpoked = Eunpoked2D[mp.Fend.corr.maskBool]
46-
DeltaV = 1e-4
47-
G1slow = np.zeros((mp.Fend.corr.Npix, Nind), dtype=complex)
48-
G2slow = np.zeros((mp.Fend.corr.Npix, Nind), dtype=complex)
49-
50-
for ii in range(Nind):
51-
# DM1
52-
mp.dm1.V = np.zeros((mp.dm1.Nact, mp.dm1.Nact))
53-
mp.dm2.V = np.zeros((mp.dm2.Nact, mp.dm2.Nact))
54-
mp.dm1.V[np.unravel_index(indG1subset[ii], mp.dm1.V.shape)] = DeltaV
55-
Epoked2D = falco.model.compact(mp, modvar)
56-
Epoked = Epoked2D[mp.Fend.corr.maskBool]
57-
G1slow[:, ii] = (Epoked - Eunpoked) / DeltaV
58-
59-
# DM2
60-
mp.dm1.V = np.zeros((mp.dm1.Nact, mp.dm1.Nact))
61-
mp.dm2.V = np.zeros((mp.dm2.Nact, mp.dm2.Nact))
62-
mp.dm2.V[np.unravel_index(indG2subset[ii], mp.dm2.V.shape)] = DeltaV
63-
Epoked2D = falco.model.compact(mp, modvar)
64-
Epoked = Epoked2D[mp.Fend.corr.maskBool]
65-
G2slow[:, ii] = (Epoked - Eunpoked) / DeltaV
6613

67-
# Tests
68-
rmsNormErrorDM1 = (np.sqrt(np.sum(np.abs(G1slow - G1fast)**2) /
69-
np.sum(np.abs(G1slow)**2)))
70-
rmsNormErrorDM2 = (np.sqrt(np.sum(np.abs(G2slow - G2fast)**2) /
71-
np.sum(np.abs(G2slow)**2)))
72-
73-
print('rmsNormErrorDM1 = %.3f' % rmsNormErrorDM1)
74-
print('rmsNormErrorDM2 = %.3f' % rmsNormErrorDM2)
75-
assert rmsNormErrorDM1 < 0.01
76-
assert rmsNormErrorDM2 < 0.01
14+
# Create a Generator instance with a seed
15+
rng1 = np.random.default_rng(seed=42)
16+
rng2 = np.random.default_rng(seed=7)
17+
18+
for case_num in range(6):
19+
20+
mp = deepcopy(CONFIG.mp)
21+
mp.runLabel = 'test_vc'
22+
23+
if case_num == 0:
24+
mp.flagRotation = False # Whether to rotate 180 degrees between conjugate planes in the compact model
25+
mp.NrelayFend = 0 # How many times to rotate the final image by 180 degrees
26+
elif case_num == 1:
27+
mp.flagRotation = True # Whether to rotate 180 degrees between conjugate planes in the compact model
28+
mp.Nrelay1to2 = 0
29+
mp.Nrelay2to3 = 0
30+
mp.Nrelay3to4 = 0
31+
mp.NrelayFend = 0 # How many times to rotate the final image by 180 degrees
32+
elif case_num == 2:
33+
mp.flagRotation = True # Whether to rotate 180 degrees between conjugate planes in the compact model
34+
mp.Nrelay1to2 = 1
35+
mp.Nrelay2to3 = 0
36+
mp.Nrelay3to4 = 0
37+
mp.NrelayFend = 0 # How many times to rotate the final image by 180 degrees
38+
elif case_num == 3:
39+
mp.flagRotation = True # Whether to rotate 180 degrees between conjugate planes in the compact model
40+
mp.Nrelay1to2 = 0
41+
mp.Nrelay2to3 = 1
42+
mp.Nrelay3to4 = 0
43+
mp.NrelayFend = 0 # How many times to rotate the final image by 180 degrees
44+
elif case_num == 4:
45+
mp.flagRotation = True # Whether to rotate 180 degrees between conjugate planes in the compact model
46+
mp.Nrelay1to2 = 0
47+
mp.Nrelay2to3 = 0
48+
mp.Nrelay3to4 = 1
49+
mp.NrelayFend = 0 # How many times to rotate the final image by 180 degrees
50+
elif case_num == 5:
51+
mp.flagRotation = True # Whether to rotate 180 degrees between conjugate planes in the compact model
52+
mp.Nrelay1to2 = 0
53+
mp.Nrelay2to3 = 0
54+
mp.Nrelay3to4 = 0
55+
mp.NrelayFend = 1 # How many times to rotate the final image by 180 degrees
56+
else:
57+
raise ValueError('Case not defined.')
58+
59+
print('\n*** NEW TEST CASE ***')
60+
print(f'mp.flagRotation = {mp.flagRotation}')
61+
print(f'mp.Nrelay1to2 = {mp.Nrelay1to2}')
62+
print(f'mp.Nrelay2to3 = {mp.Nrelay2to3}')
63+
print(f'mp.Nrelay3to4 = {mp.Nrelay3to4}')
64+
print(f'mp.NrelayFend = {mp.NrelayFend}')
65+
print('')
66+
67+
mp = deepcopy(CONFIG.mp)
68+
mp.runLabel = 'test_lc'
69+
70+
scalefac = 1 # 0.5
71+
dm1v0 = scalefac*(rng1.random((mp.dm1.Nact, mp.dm1.Nact))-0.5)
72+
dm2v0 = scalefac*(rng2.random((mp.dm2.Nact, mp.dm2.Nact))-0.5)
73+
74+
mp.path = falco.config.Object()
75+
LOCAL_PATH = os.path.dirname(os.path.abspath(__file__))
76+
mp.path.falco = os.path.dirname(os.path.dirname(LOCAL_PATH))
77+
_ = falco.setup.flesh_out_workspace(mp)
78+
79+
# Fast Jacobian calculation
80+
mp.dm1.V = dm1v0.copy() # np.zeros((mp.dm1.Nact, mp.dm1.Nact))
81+
mp.dm2.V = dm2v0.copy() # np.zeros((mp.dm2.Nact, mp.dm2.Nact))
82+
jacStruct = falco.model.jacobian(mp) # Get structure containing Jacobians
83+
84+
G1fastAll = jacStruct.G1
85+
G2fastAll = jacStruct.G2
86+
Nind = 20
87+
subinds = np.ix_(np.arange(3, 20*4, 4).astype(int))
88+
absG1sum = np.sum(np.abs(G1fastAll), axis=1)
89+
indG1 = np.nonzero(absG1sum > 1e-2*np.max(absG1sum))[0]
90+
indG1subset = indG1[subinds] # Take a 20-actuator subset
91+
absG2sum = np.sum(np.abs(G2fastAll), axis=1)
92+
indG2 = np.nonzero(absG2sum > 1e-2*np.max(absG2sum))[0]
93+
indG2subset = indG2[subinds] # Take a 20-actuator subset
94+
G1fast = np.squeeze(G1fastAll[:, indG1subset])
95+
G2fast = np.squeeze(G2fastAll[:, indG2subset])
96+
97+
# Compute Jacobian via differencing (slower)
98+
modvar = falco.config.ModelVariables()
99+
modvar.whichSource = 'star'
100+
modvar.sbpIndex = 0
101+
modvar.starIndex = 0
102+
Eunpoked2D = falco.model.compact(mp, modvar)
103+
Eunpoked = Eunpoked2D[mp.Fend.corr.maskBool]
104+
DeltaV = 1e-4
105+
G1slow = np.zeros((mp.Fend.corr.Npix, Nind), dtype=complex)
106+
G2slow = np.zeros((mp.Fend.corr.Npix, Nind), dtype=complex)
107+
108+
for ii in range(Nind):
109+
# DM1
110+
mp.dm1.V = dm1v0.copy() # np.zeros((mp.dm1.Nact, mp.dm1.Nact))
111+
mp.dm2.V = dm2v0.copy() # np.zeros((mp.dm2.Nact, mp.dm2.Nact))
112+
mp.dm1.V[np.unravel_index(indG1subset[ii], mp.dm1.V.shape)] += DeltaV
113+
Epoked2D = falco.model.compact(mp, modvar)
114+
Epoked = Epoked2D[mp.Fend.corr.maskBool]
115+
G1slow[:, ii] = (Epoked - Eunpoked) / DeltaV
116+
117+
# DM2
118+
mp.dm1.V = dm1v0.copy() # np.zeros((mp.dm1.Nact, mp.dm1.Nact))
119+
mp.dm2.V = dm2v0.copy() # np.zeros((mp.dm2.Nact, mp.dm2.Nact))
120+
mp.dm2.V[np.unravel_index(indG2subset[ii], mp.dm2.V.shape)] += DeltaV
121+
Epoked2D = falco.model.compact(mp, modvar)
122+
Epoked = Epoked2D[mp.Fend.corr.maskBool]
123+
G2slow[:, ii] = (Epoked - Eunpoked) / DeltaV
124+
125+
# Tests
126+
rmsNormErrorDM1 = (np.sqrt(np.sum(np.abs(G1slow - G1fast)**2) /
127+
np.sum(np.abs(G1slow)**2)))
128+
rmsNormErrorDM2 = (np.sqrt(np.sum(np.abs(G2slow - G2fast)**2) /
129+
np.sum(np.abs(G2slow)**2)))
130+
131+
print('rmsNormErrorDM1 = %.3f' % rmsNormErrorDM1)
132+
print('rmsNormErrorDM2 = %.3f' % rmsNormErrorDM2)
133+
assert rmsNormErrorDM1 < 0.01
134+
assert rmsNormErrorDM2 < 0.01
77135

78136

79137
def test_jacobian_lc_no_fpm():

‎tests/functional/test_jacobian_vc.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -73,7 +73,7 @@ def test_jacobian_vc_fft():
7373
mp.jac.mftToVortex = False
7474

7575
scalefac = 1 # 0.5
76-
dm1v0 = scalefac*(rng1.random((mp.dm1.Nact, mp.dm1.Nact))-0.5)
76+
dm1v0 = scalefac*(rng1.random((mp.dm1.Nact, mp.dm1.Nact))-0.5)
7777
dm2v0 = scalefac*(rng2.random((mp.dm2.Nact, mp.dm2.Nact))-0.5)
7878

7979
inputs = {

0 commit comments

Comments
 (0)