From d73ca5852ce7622c59c9dbfb6a7a196dd9837e5f Mon Sep 17 00:00:00 2001 From: A J Eldorado Riggs Date: Mon, 10 Nov 2025 10:48:31 -0800 Subject: [PATCH 1/3] Add spline resample and testing scripts --- falco/diff_dm.py | 128 ++++++++- falco/est.py | 6 +- falco/util.py | 2 +- tests/functional/check_gradient_both_dms.py | 5 +- tests/functional/check_probe_amplitude.py | 287 ++++++++++++++++++++ tests/functional/check_resampling.py | 287 ++++++++++++++++++++ tests/functional/test_ad_efc_gradient_lc.py | 4 +- tests/unit/test_offcenter_crop.py | 198 ++++++++++++++ 8 files changed, 906 insertions(+), 11 deletions(-) create mode 100644 tests/functional/check_probe_amplitude.py create mode 100644 tests/functional/check_resampling.py create mode 100644 tests/unit/test_offcenter_crop.py diff --git a/falco/diff_dm.py b/falco/diff_dm.py index 4290949..80367bd 100644 --- a/falco/diff_dm.py +++ b/falco/diff_dm.py @@ -4,11 +4,18 @@ import os +from scipy.interpolate import RectBivariateSpline from scipy import ndimage +try: + from scipy.ndimage import map_coordinates +except: + from scipy.ndimage.interpolation import map_coordinates -from falco import fftutils, util, proper +from falco import fftutils, util, proper, mask, check from astropy.io import fits + + # the functions: # make_rotation_matrix # promote_3d_transformation_to_homography @@ -417,7 +424,120 @@ def prepare_actuator_lattice(shape, Nact, sep, dtype): } +def spline_resample(array_in, zoom): + """ + Translate, rotate, and downsample a pixel-centered mask. + + Parameters + ---------- + array_in : np.ndarray + 2-D, pixel-centered array containing the pupil mask. + zoom : int/float or tuple of two ints/floats + + Returns + ------- + arrayOut : np.ndarray + 2-D, even-sized, square array containing the resized mask. + """ + check.twoD_array(array_in, 'array_in', TypeError) + # check.real_positive_scalar(nBeamIn, 'nBeamIn', TypeError) + + # if mag <= 0: + # raise ValueError('This function is for downsampling only.') + # else: + # if array_in.shape[0] % 2 == 0: # if in an even-sized array + # # Crop assuming mask is pixel centered with row 0 and column 0 empty. + # array_in = array_in[1::, 1::] + + if zoom == 1: + return array_in + + if isinstance(zoom, (float, int)): + zoom = (zoom, zoom) + elif not isinstance(zoom, tuple): + zoom = tuple(float(zoom) for zoom in zoom) + + array_in = util.pad_crop(array_in, int(np.max(array_in.shape))) + + # Array sizes + narray_in = array_in.shape[0] + narray_out = util.ceil_odd(narray_in*np.max(zoom)+2) + # 2 pixels added to guarantee the offset mask is fully contained in + # the output array. + nBeamIn = 1 # Only relative values matter + nBeamOutX = nBeamIn*zoom[1] + nBeamOutY = nBeamIn*zoom[0] + dxIn = 1./nBeamIn + dYin = 1./nBeamIn + dxOut = 1./nBeamOutX + dyOut = 1./nBeamOutY + + # array-centered coordinates of input matrix [pupil diameters] + x0 = np.arange(-(narray_in-1.)/2., (narray_in)/2., 1)*dxIn + y0 = np.arange(-(narray_in-1.)/2., (narray_in)/2., 1)*dYin + + if False: #zoom[0] < 1 and zoom[1] < 1: + [X0, Y0] = np.meshgrid(x0, y0) + R2 = X0**2 + Y0**2 + Window = np.zeros_like(R2) + Window[R2 <= (dxOut/2.)**2 + (dyOut/2.)**2] = 1 + # R0 = np.sqrt(X0**2 + Y0**2) + # Window = 0*R0 + # Window[R0 <= dxOut/2.] = 1 + Window = Window/np.sum(Window) + + # To get good grayscale edges, convolve with the correct window + # before downsampling. + f_window = np.fft.ifft2(np.fft.ifftshift(Window))*narray_in + f_array_in = np.fft.ifft2(np.fft.ifftshift(array_in))*narray_in + A = np.fft.fftshift(np.fft.fft2(f_window*f_array_in)) + A = np.real(A) + else: + A = array_in + + x1 = (np.arange(-(narray_out-1.)/2., narray_out/2., 1) - 0)*dxOut + y1 = (np.arange(-(narray_out-1.)/2., narray_out/2., 1) - 0)*dyOut + + # RectBivariateSpline is faster in 2-D than interp2d + interp_spline = RectBivariateSpline(y0, x0, A) + Atemp = interp_spline(y1, x1) + + arrayOut = Atemp + + # arrayOut = np.zeros((narray_out+1, narray_out+1)) + # arrayOut[1::, 1::] = np.real(Atemp) + + return arrayOut + + +def map_resample(f, zoom): + print(f'zoom = {zoom}') + if zoom == 1: + return f + + if isinstance(zoom, (float, int)): + zoom = (zoom, zoom) + elif not isinstance(zoom, tuple): + zoom = tuple(float(zoom) for zoom in zoom) + + N0 = util.ceil_odd(np.max(f.shape) + np.max(zoom) + 1) + f = util.pad_crop(f, N0) + + M = int(util.ceil_odd(N0*zoom[0]+1)) + N = int(util.ceil_odd(N0*zoom[1]+1)) + + x = (np.arange(N, dtype=np.float64) - np.floor(N/2))/zoom[1] + np.floor(N0/2) + y = (np.arange(M, dtype=np.float64) - np.floor(M/2))/zoom[0] + np.floor(N0/2) + xygrid = np.meshgrid(x, y) + fprime = map_coordinates(f.T, xygrid, order=3, mode="nearest") + + # fprime = proper.prop_magnify(f, zoom, size_out0=N1, AMP_CONSERVE=True) + + return fprime + + def fourier_resample(f, zoom): + print(f'zoom = {zoom}') if zoom == 1: return f @@ -676,7 +796,8 @@ def render(self, Nout=None, wfe=True): warped *= (2*self.obliquity) if self.upsample != 1: - warped = fourier_resample(warped, self.upsample) + # warped = fourier_resample(warped, self.upsample) + warped = spline_resample(warped, self.upsample) self.Nintermediate = warped.shape warped = util.pad_crop(warped, Nout) @@ -724,7 +845,8 @@ def render_backprop(self, protograd, wfe=True): protograd = util.pad_crop(protograd, self.Nintermediate) if self.upsample != 1: upsample = self.ifn.shape[0]/protograd.shape[0] - protograd = fourier_resample(protograd, upsample) + # protograd = fourier_resample(protograd, upsample) + protograd = spline_resample(protograd, upsample) protograd /= upsample**2 if wfe: diff --git a/falco/est.py b/falco/est.py index 6e38f8e..c0a956e 100644 --- a/falco/est.py +++ b/falco/est.py @@ -367,7 +367,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])): for iStar in range(mp.compact.star.count): - modvar = falco.config.ModelVariables() # Initialize the new structure + modvar = falco.config.ModelVariables() # Initialize the new structure modvar.starIndex = iStar modvar.whichSource = 'star' @@ -375,7 +375,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])): modvar.sbpIndex = iSubband print('Wavelength: %u/%u ... ' % (iSubband, mp.Nsbp-1)) - modeIndex = iStar*mp.Nsbp + iSubband #(iStar-1)*mp.Nsbp + iSubband + modeIndex = iStar*mp.Nsbp + iSubband print('Mode: %u/%u ... ' % (modeIndex, mp.jac.Nmode-1)) # Measure current contrast level average @@ -503,7 +503,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])): ampSq = (Iplus+Iminus)/2 - np.tile(I0vec.reshape((-1, 1)), (1, Npairs)) # square of probe E-field amplitudes 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) + # isnonzero = np.all(amp, 1) 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 diff --git a/falco/util.py b/falco/util.py index 00fd321..89be154 100644 --- a/falco/util.py +++ b/falco/util.py @@ -521,7 +521,7 @@ def offcenter_crop(arrayIn, centerRow, centerCol, nRowOut, nColOut): x_pad = int(np.max((colPadPre, colPadPost))) arrayPadded = pad_crop(arrayIn, (nRowIn+2*y_pad, - nColIn+2*x_pad)) + nColIn+2*x_pad)) centerCol += x_pad centerRow += y_pad diff --git a/tests/functional/check_gradient_both_dms.py b/tests/functional/check_gradient_both_dms.py index 0b9d650..3b9f54c 100644 --- a/tests/functional/check_gradient_both_dms.py +++ b/tests/functional/check_gradient_both_dms.py @@ -10,7 +10,7 @@ from falco.util import pad_crop -show_plots = False +show_plots = True n_pupil = 50 n_pupil_total = n_pupil*n_pupil @@ -108,7 +108,7 @@ def jac_model(masks, u1, u2, which_dm, du): elif which_dm == 2: e_dm2 = falco.prop.ptp(e_dm1, dx_pupil_m*n_pupil_pad, lam, d_dm1_dm2) e_dm2 *= np.exp(1j*kk*u2_pad) - de_dm2 = (1j*kk*du_pad) * e_dm2 + de_dm2 = (1j*kk*du_pad) * e_dm2 de_dm1_eff = falco.prop.ptp(de_dm2, dx_pupil_m*n_pupil_pad, lam, -d_dm1_dm2) else: @@ -158,7 +158,6 @@ def reverse_model(masks, u1, u2, du1, du2, Eest2D, I00): e_dm2_grad *= np.conj(np.exp(1j*kk*u2_pad)) phase_dm2_bar = -kk*np.imag(e_dm2_grad * np.conj(e_dm2_pre)) - e_dm1_grad = falco.prop.ptp(e_dm2_grad, dx_pupil_m*n_pupil_pad, lam, -d_dm1_dm2) phase_dm1_bar = -kk*np.imag(e_dm1_grad * np.conj(e_dm1_post)) diff --git a/tests/functional/check_probe_amplitude.py b/tests/functional/check_probe_amplitude.py new file mode 100644 index 0000000..a6e97c6 --- /dev/null +++ b/tests/functional/check_probe_amplitude.py @@ -0,0 +1,287 @@ +"""Jacobian accuracy tests for the Lyot coronagraph.""" +from copy import deepcopy +import os +from types import SimpleNamespace + +import matplotlib.pyplot as plt +import numpy as np + +import falco +from falco.util import pad_crop, ceil_odd +from falco.diff_dm import spline_resample, map_resample, fourier_resample + + +show_plots = True + +n_pupil = 50 +n_pupil_total = n_pupil*n_pupil +n_pupil_pad = 100 +centering = 'pixel' + +res_cam = 3 #3 +rho0 = 2.5 +rho1 = 9.0 +n_cam = falco.util.ceil_even(2*rho1*res_cam + 1) + +fl = 1 +lam_norm = 1 +diam_pupil_m = 46e-3 +dx_pupil = 1/n_pupil +dxi_cam = 1/res_cam + +lam = 575e-9 # meters +gain = 1e-9 +kk = gain*2*np.pi/lam + + + + +def make_masks(): + + # masks = SimpleNamespace() #falco.config.Object() + masks = falco.config.Object() + + # Make pupil mask + + inputs = { + 'Nbeam': n_pupil, + 'Npad': n_pupil_total, + 'OD': 1, + 'wStrut': 0.04, + # 'centering': 'interpixel', + # 'angStrut': 10 + np.array([0, 120, 240]) + } + pupil = falco.mask.falco_gen_pupil_Simple(inputs) + + + # pupil = falco.mask.falco_gen_pupil_Roman_CGI_20200513(n_pupil, centering) + pupil = falco.util.pad_crop(pupil, n_pupil_pad) + masks.pupil = pupil + + # Make software mask at detector + SCORE = {} + SCORE["Nxi"] = n_cam + SCORE["Neta"] = n_cam + SCORE["pixresFP"] = res_cam + SCORE["centering"] = centering + SCORE["rhoInner"] = rho0 # lambda0/D + SCORE["rhoOuter"] = rho1 # lambda0/D + SCORE["angDeg"] = 180 # degrees + SCORE["whichSide"] = 'both' + SCORE["shape"] = 'circle' + [maskScore, _, _] = falco.mask.falco_gen_SW_mask(SCORE) + masks.dh = np.array(maskScore, dtype=bool) + + return masks + + +def make_jacobian(masks, u1): + + n_pix = np.sum(masks.dh.astype(int)) + jac = np.zeros((n_pix, n_pupil_total), dtype=complex) + + du1_flat = np.zeros((n_pupil_total,)) + for ii in range(n_pupil_total): + du1_flat[ii] = 1 + du1 = du1_flat.reshape((n_pupil, n_pupil)) + + jac[:, ii] = jac_model(masks, u1, du1) + + du1_flat *= 0 # reset + + return jac + + +def jac_model(masks, u1, du1): + + u1_pad = pad_crop(u1, n_pupil_pad) + du1_pad = pad_crop(du1, n_pupil_pad) + de_pupil = masks.pupil * np.exp(1j*kk*u1_pad) * (1j*kk*du1_pad) + + de_cam = falco.prop.mft_p2f(de_pupil, fl, lam_norm, dx_pupil, dxi_cam, n_cam, dxi_cam, n_cam, centering) + de_vec = de_cam[masks.dh] + + return de_vec + + +def forward_model(masks, u1): + + u1_pad = falco.util.pad_crop(u1, n_pupil_pad) + e_dm1_post = np.exp(1j*kk*u1_pad) * masks.pupil + + e_cam = falco.prop.mft_p2f(e_dm1_post, fl, lam_norm, dx_pupil, dxi_cam, n_cam, dxi_cam, n_cam, centering) + e_vec = e_cam[masks.dh] + + return e_vec, e_cam + + +def reverse_model(masks, u1, du1, Eest2D, I00): + + _, EFendA = forward_model(masks, u1) + _, EFendB = forward_model(masks, u1+du1) + dEend = EFendB - EFendA + EdhNew = Eest2D + dEend/np.sqrt(I00) + + u1_pad = falco.util.pad_crop(u1, n_pupil_pad) + e_dm1_post = np.exp(1j*kk*u1_pad) * masks.pupil + # e_dm1_post = masks.pupil # DEBUGGING + + # Gradient + # cam_grad = EdhNew*np.real(masks.dh.astype(float)) / np.sum(np.abs(Eest2D[masks.dh])**2) + cam_grad = 2/np.sqrt(I00)*EdhNew*masks.dh #np.real(masks.dh.astype(float)) + # cam_grad = np.zeros_like(EdhNew) + # cam_grad[masks.dh] = 2*EdhNew[masks.dh] + + pupil_grad = falco.prop.mft_f2p(cam_grad, -fl, lam_norm, dxi_cam, dxi_cam, dx_pupil, n_pupil_pad, centering) + + phase_DM2_bar = -kk*np.imag(pupil_grad * np.conj(e_dm1_post)) + gradient = pad_crop(phase_DM2_bar, n_pupil).ravel() + + return gradient + + +def test_resampling(): + masks = make_masks() + + usfac = 2.9 + dsfac = 0.33 + + pupil0 = masks.pupil + pupil0 = pad_crop(pupil0, int(ceil_odd(pupil0.shape[0]))) + + pupil_spline_up = spline_resample(pupil0, usfac) + pupil_spline_down = spline_resample(pupil0, dsfac) + # pupil_spline_up = map_resample(pupil0, usfac) + # pupil_spline_down = map_resample(pupil0, dsfac) + + pupil_fourier_up = fourier_resample(pupil0, usfac) + pupil_fourier_down = fourier_resample(pupil0, dsfac) + + pupil_fourier_up = pad_crop(pupil_fourier_up, pupil_spline_up.shape) + pupil_fourier_down = pad_crop(pupil_fourier_down, pupil_spline_down.shape) + + print(f'pupil0.shape = {pupil0.shape}') + print(f'pupil_spline_up.shape = {pupil_spline_up.shape}') + print(f'pupil_spline_down.shape = {pupil_spline_down.shape}') + # print(f'pupil_fourier_up.shape = {pupil_fourier_up.shape}') + # print(f'pupil_fourier_down.shape = {pupil_fourier_down.shape}') + + if show_plots: + + plt.figure(1) + plt.imshow(pupil0) + plt.colorbar() + plt.title('pupil0') + + plt.figure(2) + plt.imshow(pupil_spline_up) + plt.colorbar() + plt.title('pupil_spline_up') + + plt.figure(3) + plt.imshow(pupil_spline_down) + plt.colorbar() + plt.title('pupil_spline_down') + + plt.figure(4) + plt.imshow(pupil_fourier_up - pupil_spline_up) + plt.colorbar() + plt.title('pupil_fourier_up') + + plt.figure(5) + plt.imshow(pupil_fourier_down - pupil_spline_down) + plt.colorbar() + plt.title('pupil_fourier_down') + + + plt.show() + + +def test_adjoint(): + """Most basic tests of the adjoint model accuracy.""" + + masks = make_masks() + # u1 = np.zeros((n_pupil, n_pupil)) + u1 = 0.8 * (np.random.rand(n_pupil, n_pupil)-0.5) + + # Make Jacobian G1 + G1 = make_jacobian(masks, u1) + + # Get true E-field, e_vec + e_vec, _ = forward_model(masks, u1) + + # Get "estimated" starting E-field + _, Eest2D = forward_model(masks, u1) + I00 = np.max(np.abs(Eest2D)**2) + # I00 = 1 # DEBUGGING + print(f'I00 = {I00}') + print(f'sqrt(I00) = {np.sqrt(I00)}') + print(f'kk = {kk}') + + e_vec /= np.sqrt(I00) + Eest2D /= np.sqrt(I00) + G1 /= np.sqrt(I00) + + # Assign delta command + du1 = np.eye(n_pupil) + + # Compute expected response + u1_bar_expected = 2 * G1.conj().T @ (e_vec + G1 @ du1.ravel()) + u1_bar_expected_2d = u1_bar_expected.reshape((n_pupil, n_pupil)) + + # Compute gradient + gradient = reverse_model(masks, u1, du1, Eest2D, I00) + + u1_bar_out_2d = gradient.reshape((n_pupil, n_pupil)) + + # breakpoint() + + if show_plots: + + plt.figure(1) + plt.imshow(np.real(u1_bar_expected_2d)) + plt.colorbar() + plt.title('u1_bar_expected_2d, real') + + plt.figure(11) + plt.imshow(np.real(u1_bar_out_2d)) + plt.colorbar() + plt.title('u1_bar_out_2d, real') + + # plt.figure(21) + # plt.imshow(np.real(u1_bar_expected_2d - u1_bar_out_2d)) + # plt.colorbar() + # plt.title('u1_bar_expected_2d - u1_bar_out_2d, real') + + plt.figure(31) + plt.imshow(np.real(u1_bar_expected_2d / u1_bar_out_2d)) + plt.colorbar() + plt.title('u1_bar_expected_2d / u1_bar_out_2d, real') + + # plt.figure(2) + # plt.imshow(np.real(u2_bar_expected_2d)) + # plt.colorbar() + # plt.title('u2_bar_expected_2d, real') + + # plt.figure(12) + # plt.imshow(np.real(u2_bar_out_2d)) + # plt.colorbar() + # plt.title('u2_bar_out_2d, real') + + + plt.figure(10) + plt.imshow(u1) + plt.colorbar() + plt.title('u1') + + plt.figure(20) + plt.imshow(np.log10(np.abs(Eest2D))) + plt.colorbar() + plt.title('log10(abs(Eest2D))') + + plt.show() + + +if __name__ == '__main__': + # test_adjoint() + test_resampling() diff --git a/tests/functional/check_resampling.py b/tests/functional/check_resampling.py new file mode 100644 index 0000000..a6e97c6 --- /dev/null +++ b/tests/functional/check_resampling.py @@ -0,0 +1,287 @@ +"""Jacobian accuracy tests for the Lyot coronagraph.""" +from copy import deepcopy +import os +from types import SimpleNamespace + +import matplotlib.pyplot as plt +import numpy as np + +import falco +from falco.util import pad_crop, ceil_odd +from falco.diff_dm import spline_resample, map_resample, fourier_resample + + +show_plots = True + +n_pupil = 50 +n_pupil_total = n_pupil*n_pupil +n_pupil_pad = 100 +centering = 'pixel' + +res_cam = 3 #3 +rho0 = 2.5 +rho1 = 9.0 +n_cam = falco.util.ceil_even(2*rho1*res_cam + 1) + +fl = 1 +lam_norm = 1 +diam_pupil_m = 46e-3 +dx_pupil = 1/n_pupil +dxi_cam = 1/res_cam + +lam = 575e-9 # meters +gain = 1e-9 +kk = gain*2*np.pi/lam + + + + +def make_masks(): + + # masks = SimpleNamespace() #falco.config.Object() + masks = falco.config.Object() + + # Make pupil mask + + inputs = { + 'Nbeam': n_pupil, + 'Npad': n_pupil_total, + 'OD': 1, + 'wStrut': 0.04, + # 'centering': 'interpixel', + # 'angStrut': 10 + np.array([0, 120, 240]) + } + pupil = falco.mask.falco_gen_pupil_Simple(inputs) + + + # pupil = falco.mask.falco_gen_pupil_Roman_CGI_20200513(n_pupil, centering) + pupil = falco.util.pad_crop(pupil, n_pupil_pad) + masks.pupil = pupil + + # Make software mask at detector + SCORE = {} + SCORE["Nxi"] = n_cam + SCORE["Neta"] = n_cam + SCORE["pixresFP"] = res_cam + SCORE["centering"] = centering + SCORE["rhoInner"] = rho0 # lambda0/D + SCORE["rhoOuter"] = rho1 # lambda0/D + SCORE["angDeg"] = 180 # degrees + SCORE["whichSide"] = 'both' + SCORE["shape"] = 'circle' + [maskScore, _, _] = falco.mask.falco_gen_SW_mask(SCORE) + masks.dh = np.array(maskScore, dtype=bool) + + return masks + + +def make_jacobian(masks, u1): + + n_pix = np.sum(masks.dh.astype(int)) + jac = np.zeros((n_pix, n_pupil_total), dtype=complex) + + du1_flat = np.zeros((n_pupil_total,)) + for ii in range(n_pupil_total): + du1_flat[ii] = 1 + du1 = du1_flat.reshape((n_pupil, n_pupil)) + + jac[:, ii] = jac_model(masks, u1, du1) + + du1_flat *= 0 # reset + + return jac + + +def jac_model(masks, u1, du1): + + u1_pad = pad_crop(u1, n_pupil_pad) + du1_pad = pad_crop(du1, n_pupil_pad) + de_pupil = masks.pupil * np.exp(1j*kk*u1_pad) * (1j*kk*du1_pad) + + de_cam = falco.prop.mft_p2f(de_pupil, fl, lam_norm, dx_pupil, dxi_cam, n_cam, dxi_cam, n_cam, centering) + de_vec = de_cam[masks.dh] + + return de_vec + + +def forward_model(masks, u1): + + u1_pad = falco.util.pad_crop(u1, n_pupil_pad) + e_dm1_post = np.exp(1j*kk*u1_pad) * masks.pupil + + e_cam = falco.prop.mft_p2f(e_dm1_post, fl, lam_norm, dx_pupil, dxi_cam, n_cam, dxi_cam, n_cam, centering) + e_vec = e_cam[masks.dh] + + return e_vec, e_cam + + +def reverse_model(masks, u1, du1, Eest2D, I00): + + _, EFendA = forward_model(masks, u1) + _, EFendB = forward_model(masks, u1+du1) + dEend = EFendB - EFendA + EdhNew = Eest2D + dEend/np.sqrt(I00) + + u1_pad = falco.util.pad_crop(u1, n_pupil_pad) + e_dm1_post = np.exp(1j*kk*u1_pad) * masks.pupil + # e_dm1_post = masks.pupil # DEBUGGING + + # Gradient + # cam_grad = EdhNew*np.real(masks.dh.astype(float)) / np.sum(np.abs(Eest2D[masks.dh])**2) + cam_grad = 2/np.sqrt(I00)*EdhNew*masks.dh #np.real(masks.dh.astype(float)) + # cam_grad = np.zeros_like(EdhNew) + # cam_grad[masks.dh] = 2*EdhNew[masks.dh] + + pupil_grad = falco.prop.mft_f2p(cam_grad, -fl, lam_norm, dxi_cam, dxi_cam, dx_pupil, n_pupil_pad, centering) + + phase_DM2_bar = -kk*np.imag(pupil_grad * np.conj(e_dm1_post)) + gradient = pad_crop(phase_DM2_bar, n_pupil).ravel() + + return gradient + + +def test_resampling(): + masks = make_masks() + + usfac = 2.9 + dsfac = 0.33 + + pupil0 = masks.pupil + pupil0 = pad_crop(pupil0, int(ceil_odd(pupil0.shape[0]))) + + pupil_spline_up = spline_resample(pupil0, usfac) + pupil_spline_down = spline_resample(pupil0, dsfac) + # pupil_spline_up = map_resample(pupil0, usfac) + # pupil_spline_down = map_resample(pupil0, dsfac) + + pupil_fourier_up = fourier_resample(pupil0, usfac) + pupil_fourier_down = fourier_resample(pupil0, dsfac) + + pupil_fourier_up = pad_crop(pupil_fourier_up, pupil_spline_up.shape) + pupil_fourier_down = pad_crop(pupil_fourier_down, pupil_spline_down.shape) + + print(f'pupil0.shape = {pupil0.shape}') + print(f'pupil_spline_up.shape = {pupil_spline_up.shape}') + print(f'pupil_spline_down.shape = {pupil_spline_down.shape}') + # print(f'pupil_fourier_up.shape = {pupil_fourier_up.shape}') + # print(f'pupil_fourier_down.shape = {pupil_fourier_down.shape}') + + if show_plots: + + plt.figure(1) + plt.imshow(pupil0) + plt.colorbar() + plt.title('pupil0') + + plt.figure(2) + plt.imshow(pupil_spline_up) + plt.colorbar() + plt.title('pupil_spline_up') + + plt.figure(3) + plt.imshow(pupil_spline_down) + plt.colorbar() + plt.title('pupil_spline_down') + + plt.figure(4) + plt.imshow(pupil_fourier_up - pupil_spline_up) + plt.colorbar() + plt.title('pupil_fourier_up') + + plt.figure(5) + plt.imshow(pupil_fourier_down - pupil_spline_down) + plt.colorbar() + plt.title('pupil_fourier_down') + + + plt.show() + + +def test_adjoint(): + """Most basic tests of the adjoint model accuracy.""" + + masks = make_masks() + # u1 = np.zeros((n_pupil, n_pupil)) + u1 = 0.8 * (np.random.rand(n_pupil, n_pupil)-0.5) + + # Make Jacobian G1 + G1 = make_jacobian(masks, u1) + + # Get true E-field, e_vec + e_vec, _ = forward_model(masks, u1) + + # Get "estimated" starting E-field + _, Eest2D = forward_model(masks, u1) + I00 = np.max(np.abs(Eest2D)**2) + # I00 = 1 # DEBUGGING + print(f'I00 = {I00}') + print(f'sqrt(I00) = {np.sqrt(I00)}') + print(f'kk = {kk}') + + e_vec /= np.sqrt(I00) + Eest2D /= np.sqrt(I00) + G1 /= np.sqrt(I00) + + # Assign delta command + du1 = np.eye(n_pupil) + + # Compute expected response + u1_bar_expected = 2 * G1.conj().T @ (e_vec + G1 @ du1.ravel()) + u1_bar_expected_2d = u1_bar_expected.reshape((n_pupil, n_pupil)) + + # Compute gradient + gradient = reverse_model(masks, u1, du1, Eest2D, I00) + + u1_bar_out_2d = gradient.reshape((n_pupil, n_pupil)) + + # breakpoint() + + if show_plots: + + plt.figure(1) + plt.imshow(np.real(u1_bar_expected_2d)) + plt.colorbar() + plt.title('u1_bar_expected_2d, real') + + plt.figure(11) + plt.imshow(np.real(u1_bar_out_2d)) + plt.colorbar() + plt.title('u1_bar_out_2d, real') + + # plt.figure(21) + # plt.imshow(np.real(u1_bar_expected_2d - u1_bar_out_2d)) + # plt.colorbar() + # plt.title('u1_bar_expected_2d - u1_bar_out_2d, real') + + plt.figure(31) + plt.imshow(np.real(u1_bar_expected_2d / u1_bar_out_2d)) + plt.colorbar() + plt.title('u1_bar_expected_2d / u1_bar_out_2d, real') + + # plt.figure(2) + # plt.imshow(np.real(u2_bar_expected_2d)) + # plt.colorbar() + # plt.title('u2_bar_expected_2d, real') + + # plt.figure(12) + # plt.imshow(np.real(u2_bar_out_2d)) + # plt.colorbar() + # plt.title('u2_bar_out_2d, real') + + + plt.figure(10) + plt.imshow(u1) + plt.colorbar() + plt.title('u1') + + plt.figure(20) + plt.imshow(np.log10(np.abs(Eest2D))) + plt.colorbar() + plt.title('log10(abs(Eest2D))') + + plt.show() + + +if __name__ == '__main__': + # test_adjoint() + test_resampling() diff --git a/tests/functional/test_ad_efc_gradient_lc.py b/tests/functional/test_ad_efc_gradient_lc.py index b903b8d..a2a76ad 100644 --- a/tests/functional/test_ad_efc_gradient_lc.py +++ b/tests/functional/test_ad_efc_gradient_lc.py @@ -9,7 +9,7 @@ import config_wfsc_lc_quick as CONFIG -show_plots = False +show_plots = True def test_adjoint_model_lc(): @@ -242,11 +242,13 @@ def test_adjoint_model_lc(): plt.figure(51) plt.imshow(ratio1*mask1) plt.colorbar() + plt.clim(0.8, 1.2) plt.title('ratio1') plt.figure(52) plt.imshow(ratio2*mask2) plt.colorbar() + plt.clim(0.8, 1.2) plt.title('ratio2') plt.show() diff --git a/tests/unit/test_offcenter_crop.py b/tests/unit/test_offcenter_crop.py new file mode 100644 index 0000000..931b20c --- /dev/null +++ b/tests/unit/test_offcenter_crop.py @@ -0,0 +1,198 @@ +# Copyright 2025, by the California Institute of Technology. +# ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer +# at the California Institute of Technology. +""" +Unit tests for pad_crop.py +""" + +import unittest + +import numpy as np + +import numpy as np + +from falco.util import pad_crop + + +class TestOffcenterCrop(unittest.TestCase): + """ + Unit test suite for pad_crop() + + This will have lots of special cases to handle odd/even sizing and + truncating vs. non-truncating + """ + + # Success tests (non-truncating) + def test_insert_all_size_odd(self): + """Check insert behavior, (o,o) --> (o,o)""" + out = pad_crop(np.ones((3, 3)), (5, 5)) + test = np.array([[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]]) + self.assertTrue((out == test).all()) + + def test_insert_all_size_even(self): + """Check insert behavior, (e,e) --> (e,e)""" + out = pad_crop(np.ones((2, 2)), (4, 4)) + test = np.array([[0, 0, 0, 0], + [0, 1, 1, 0], + [0, 1, 1, 0], + [0, 0, 0, 0]]) + self.assertTrue((out == test).all()) + + def test_insert_first_index_odd_size_even(self): + """Check insert behavior, (o,e) --> (e,e)""" + out = pad_crop(np.ones((3, 2)), (4, 4)) + test = np.array([[0, 0, 0, 0], + [0, 1, 1, 0], + [0, 1, 1, 0], + [0, 1, 1, 0]]) + self.assertTrue((out == test).all()) + + def test_insert_second_index_odd_size_even(self): + """Check insert behavior, (e,o) --> (e,e)""" + out = pad_crop(np.ones((2, 3)), (4, 4)) + test = np.array([[0, 0, 0, 0], + [0, 1, 1, 1], + [0, 1, 1, 1], + [0, 0, 0, 0]]) + self.assertTrue((out == test).all()) + + def test_insert_first_index_even_size_odd(self): + """Check insert behavior, (e,o) --> (o,o)""" + out = pad_crop(np.ones((4, 3)), (5, 5)) + test = np.array([[0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]]) + self.assertTrue((out == test).all()) + + def test_insert_second_index_even_size_odd(self): + """Check insert behavior, (o,e) --> (o,o)""" + out = pad_crop(np.ones((3, 4)), (5, 5)) + test = np.array([[0, 0, 0, 0, 0], + [1, 1, 1, 1, 0], + [1, 1, 1, 1, 0], + [1, 1, 1, 1, 0], + [0, 0, 0, 0, 0]]) + self.assertTrue((out == test).all()) + + # Truncating + def test_insert_all_size_odd_trunc(self): + """Check insert behavior, (o,o) --> (o,o), truncating""" + tmat = np.outer(np.arange(0, 10, 2), np.arange(0, 15, 3)) + out = pad_crop(tmat, (3, 3)) + test = np.array([[6, 12, 18], + [12, 24, 36], + [18, 36, 54]]) + self.assertTrue((out == test).all()) + + def test_insert_all_size_even_trunc(self): + """Check insert behavior, (e,e) --> (e,e), truncating""" + tmat = np.outer(np.arange(0, 8, 2), np.arange(0, 12, 3)) + out = pad_crop(tmat, (2, 2)) + test = np.array([[6, 12], + [12, 24]]) + self.assertTrue((out == test).all()) + + def test_insert_first_index_odd_size_even_trunc(self): + """Check insert behavior, (e,e) --> (o,e), truncating""" + tmat = np.outer(np.arange(0, 8, 2), np.arange(0, 12, 3)) + out = pad_crop(tmat, (3, 2)) + test = np.array([[6, 12], + [12, 24], + [18, 36]]) + self.assertTrue((out == test).all()) + + def test_insert_second_index_odd_size_even_trunc(self): + """Check insert behavior, (e,e) --> (e,o), truncating""" + tmat = np.outer(np.arange(0, 8, 2), np.arange(0, 12, 3)) + out = pad_crop(tmat, (2, 3)) + test = np.array([[6, 12, 18], + [12, 24, 36]]) + self.assertTrue((out == test).all()) + + def test_insert_first_index_even_size_odd_trunc(self): + """Check insert behavior, (e,o) --> (o,o), truncating""" + tmat = np.outer(np.arange(0, 10, 2), np.arange(0, 15, 3)) + out = pad_crop(tmat, (4, 3)) + test = np.array([[0, 0, 0], + [6, 12, 18], + [12, 24, 36], + [18, 36, 54]]) + self.assertTrue((out == test).all()) + + def test_insert_second_index_even_size_odd_trunc(self): + """Check insert behavior, (o,e) --> (o,o), truncating""" + tmat = np.outer(np.arange(0, 10, 2), np.arange(0, 15, 3)) + out = pad_crop(tmat, (3, 4)) + test = np.array([[0, 6, 12, 18], + [0, 12, 24, 36], + [0, 18, 36, 54]]) + self.assertTrue((out == test).all()) + + # Mixed + def test_insert_first_large_second_small(self): + """Check insert behavior, truncating second axis only""" + out = pad_crop(np.ones((4, 4)), (6, 2)) + test = np.array([[0, 0], + [1, 1], + [1, 1], + [1, 1], + [1, 1], + [0, 0]]) + self.assertTrue((out == test).all()) + + def test_insert_first_small_second_large(self): + """Check insert behavior, truncating first axis only""" + out = pad_crop(np.ones((4, 4)), (2, 6)) + test = np.array([[0, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 0]]) + self.assertTrue((out == test).all()) + + # Other success + def test_transitivity(self): + """ + Verify that successive pad_crops always end up at the same location. + Should be no path dependence + """ + for midsize in [(3, 3), (3, 4), (8, 7), (6, 6), (2, 2), (8, 8)]: + out1a = pad_crop(np.ones((2, 2)), midsize) + out1b = pad_crop(out1a, (8, 8)) + out2 = pad_crop(np.ones((2, 2)), (8, 8)) + self.assertTrue((out1b == out2).all()) + + # def test_dtype_passed(self): + # """Check pad_crop maintains data type""" + # for dt in [bool, np.int32, np.float32, np.float64, + # np.complex64, np.complex128]: + # inm = np.ones((2, 2), dtype=dt) + # out = pad_crop(inm, (4, 4)) + # self.assertTrue(out.dtype == inm.dtype) + + # Failure tests + def test_arr0_2darray(self): + """Check input array type valid""" + for arr0 in [(2, 2), np.zeros((2,)), np.zeros((2, 2, 2))]: + with self.assertRaises(TypeError): + pad_crop(arr0, (4, 4)) + + def test_outsize_1or2elem_list(self): + """Check outsize is 2-element list""" + for outsize in [(4, 4, 4), [], None]: + with self.assertRaises(TypeError): + pad_crop(np.zeros((2, 2)), outsize) + + def test_outsize_has_non_positive_int(self): + """Check outsize elements are positive integers""" + for outsize in [(0, 4), (-5, 4), (6.3, 4), (4.0, 4)]: + with self.assertRaises(TypeError): + pad_crop(np.zeros((2, 2)), outsize) + + +if __name__ == '__main__': + unittest.main() From 947c8380d4303a593c0b6ebb88099fffc46a5e82 Mon Sep 17 00:00:00 2001 From: A J Eldorado Riggs Date: Mon, 10 Nov 2025 11:39:24 -0800 Subject: [PATCH 2/3] Revert changes made while adding an incomplete new feature. --- falco/diff_dm.py | 10 +++++----- tests/functional/test_ad_efc_gradient_lc.py | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/falco/diff_dm.py b/falco/diff_dm.py index 80367bd..c7f4527 100644 --- a/falco/diff_dm.py +++ b/falco/diff_dm.py @@ -511,7 +511,7 @@ def spline_resample(array_in, zoom): def map_resample(f, zoom): - print(f'zoom = {zoom}') + # print(f'zoom = {zoom}') if zoom == 1: return f @@ -796,8 +796,8 @@ def render(self, Nout=None, wfe=True): warped *= (2*self.obliquity) if self.upsample != 1: - # warped = fourier_resample(warped, self.upsample) - warped = spline_resample(warped, self.upsample) + warped = fourier_resample(warped, self.upsample) + # warped = spline_resample(warped, self.upsample) self.Nintermediate = warped.shape warped = util.pad_crop(warped, Nout) @@ -845,8 +845,8 @@ def render_backprop(self, protograd, wfe=True): protograd = util.pad_crop(protograd, self.Nintermediate) if self.upsample != 1: upsample = self.ifn.shape[0]/protograd.shape[0] - # protograd = fourier_resample(protograd, upsample) - protograd = spline_resample(protograd, upsample) + protograd = fourier_resample(protograd, upsample) + # protograd = spline_resample(protograd, upsample) protograd /= upsample**2 if wfe: diff --git a/tests/functional/test_ad_efc_gradient_lc.py b/tests/functional/test_ad_efc_gradient_lc.py index a2a76ad..2f0f60c 100644 --- a/tests/functional/test_ad_efc_gradient_lc.py +++ b/tests/functional/test_ad_efc_gradient_lc.py @@ -9,7 +9,7 @@ import config_wfsc_lc_quick as CONFIG -show_plots = True +show_plots = False def test_adjoint_model_lc(): From 7559a51835e6075d665c098aae62e8c952873a77 Mon Sep 17 00:00:00 2001 From: A J Eldorado Riggs Date: Tue, 25 Nov 2025 16:52:32 -0800 Subject: [PATCH 3/3] CI: Add Python 3.14 to testing. --- .github/workflows/pyfalco_test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/pyfalco_test.yml b/.github/workflows/pyfalco_test.yml index c31b0e9..34ce4ae 100644 --- a/.github/workflows/pyfalco_test.yml +++ b/.github/workflows/pyfalco_test.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ['3.12', ] + python-version: ['3.12', '3.14'] steps: # This step checks out a copy of your repository.