Add spline resample and testing scripts

.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.

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
 
@@ -677,6 +797,7 @@ def render(self, Nout=None, wfe=True):
 
         if self.upsample != 1:
             warped = fourier_resample(warped, self.upsample)
+            # warped = spline_resample(warped, self.upsample)
 
         self.Nintermediate = warped.shape
         warped = util.pad_crop(warped, Nout)
@@ -725,6 +846,7 @@ def render_backprop(self, protograd, wfe=True):
         if self.upsample != 1:
             upsample = self.ifn.shape[0]/protograd.shape[0]
             protograd = fourier_resample(protograd, upsample)
+            # protograd = spline_resample(protograd, upsample)
             protograd /= upsample**2
 
         if wfe:

falco/est.py

@@ -367,15 +367,15 @@ 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'
 
         for iSubband in range(mp.Nsbp):
 
             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

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

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))