AD-EFC gradient debugged

examples/WFIRST_LC_AD_EFC/EXAMPLE_config_WFIRST_LC_AD_EFC.py

@@ -148,6 +148,14 @@
 ###     ];
 ### [mp.Nitr, mp.relinItrVec, mp.gridSearchItrVec, mp.ctrl.log10regSchedIn, mp.dm_ind_sched] = falco_ctrl_EFC_schedule_generator(mp.ctrl.sched_mat);
 
+mp.ctrl.sched_mat = np.array([
+    [1, -2, 12, 1, 0],
+    [1, -2, 12, 1, 0],
+    [1, -4, 12, 1, 0],
+    [1, -4, 12, 1, 0],
+    [1, -2, 12, 1, 0],
+    ])
+mp.Nitr, mp.relinItrVec, mp.gridSearchItrVec, mp.ctrl.log10regSchedIn, mp.dm_ind_sched = falco.ctrl.efc_schedule_generator(mp.ctrl.sched_mat)
 
 ### Deformable Mirrors: Influence Functions
 ##--Influence Function Options:

examples/WFIRST_LC_AD_EFC/EXAMPLE_main_WFIRST_LC_AD_EFC.py

@@ -41,6 +41,15 @@
 # mp.dm1.surfFitMethod = 'lsq'
 # mp.dm2.surfFitMethod = 'lsq'
 
+mp.ctrl.sched_mat = np.array([
+    [1, -2, 12, 1, 0],
+    [1, -2, 12, 1, 0],
+    [1, -4, 12, 1, 0],
+    [1, -4, 12, 1, 0],
+    [1, -2, 12, 1, 0],
+    ])
+mp.Nitr, mp.relinItrVec, mp.gridSearchItrVec, mp.ctrl.log10regSchedIn, mp.dm_ind_sched = falco.ctrl.efc_schedule_generator(mp.ctrl.sched_mat)
+
 
 # %% Set up the workspace
 

examples/try_running_falco/EXAMPLE_try_running_FALCO_AD_EFC.py

@@ -33,7 +33,17 @@
 mp.Nsbp = 1            # Number of sub-bandpasses to divide the whole bandpass into for estimation and control
 mp.Nwpsbp = 1          # Number of wavelengths to used to approximate an image in each sub-bandpass
 
-mp.Nitr = 3  # Number of wavefront control iterations
+mp.dm1.useDifferentiableModel = True
+mp.dm2.useDifferentiableModel = True
+
+# mp.Nitr = 3  # Number of wavefront control iterations
+mp.ctrl.sched_mat = np.array([
+    [1, -4, 12, 1, 0],
+    [1, -4, 12, 1, 0],
+    [1, -4, 12, 1, 0],
+    ])
+mp.Nitr, mp.relinItrVec, mp.gridSearchItrVec, mp.ctrl.log10regSchedIn, mp.dm_ind_sched = falco.ctrl.efc_schedule_generator(mp.ctrl.sched_mat)
+
 
 mp.controller = 'AD-EFC'
 mp.ctrl.ad = falco.config.Object()
@@ -45,11 +55,6 @@
 
 mp.ctrl.log10regVec = np.array([-6, ])
 
-# # Use least-squares surface fitting instead of back-propagation model.
-# mp.dm1.useDifferentiableModel = False
-# mp.dm2.useDifferentiableModel = False
-# mp.dm1.surfFitMethod = 'lsq'
-# mp.dm2.surfFitMethod = 'lsq'
 
 
 # %% Perform the Wavefront Sensing and Control

falco/diff_dm.py

@@ -138,7 +138,6 @@ def dm_init_falco_wrapper(dm,dx,Narray,dm_z0, dm_xc, dm_yc, spacing=0.,**kwargs)
 
     dx_inf = dx_inf * dx_dm / dx_dm_inf  # Influence function sampling scaled
                                          # to specified DM actuator spacing
-
     #else:
     dm_z_commanded = dm_z
     s = dm_z.shape
@@ -440,6 +439,7 @@ def fourier_resample(f, zoom):
     fprime = imft2_core(F, Mx, My).real
     fprime *= np.sqrt((zoom[0]*zoom[1]))
     # fprime *= np.sqrt((zoom[0]*zoom[1]))/(np.sqrt(f.size))
+    # fprime *= np.sum(fprime)/np.sum(f)
 
     return fprime
 
@@ -682,7 +682,7 @@ def render(self, Nout=None, wfe=True):
         warped = util.pad_crop(warped, Nout)
         return warped
 
-    def render_backprop(self, protograd, gain_map, wfe=True):
+    def render_backprop(self, protograd, wfe=True):
         """Gradient backpropagation for render().
 
         Parameters
@@ -725,6 +725,7 @@ def render_backprop(self, protograd, gain_map, wfe=True):
         if self.upsample != 1:
             upsample = self.ifn.shape[0]/protograd.shape[0]
             protograd = fourier_resample(protograd, upsample)
+            protograd /= upsample**2
 
         if wfe:
             protograd *= (2*self.obliquity)
@@ -734,5 +735,6 @@ def render_backprop(self, protograd, gain_map, wfe=True):
             protograd = warp(protograd, self.invprojx, self.invprojy)
 
         # return protograd
-        in_actuator_space = apply_precomputed_transfer_function( protograd, np.conj(self.tf) )
-        return in_actuator_space[self.iyy, self.ixx] / gain_map / np.sum(self.ifn)
+        in_actuator_space = apply_precomputed_transfer_function(protograd, np.conj(self.tf))
+
+        return in_actuator_space[self.iyy, self.ixx]

falco/model/__init__.py

@@ -1 +1 @@
-from .models import *
+from .models import *

falco/setup.py

@@ -65,8 +65,9 @@ def flesh_out_workspace(mp):
     out = init_storage_arrays(mp)
 
     print('\nBeginning Trial %d of Series %d.\n' % (mp.TrialNum, mp.SeriesNum))
-    print('DM 1-to-2 Fresnel number (using radius) = ' +
-          str((mp.P2.D/2)**2/(mp.d_dm1_dm2*mp.lambda0)))
+    if mp.d_dm1_dm2*mp.lambda0 != 0:
+        print('DM 1-to-2 Fresnel number (using radius) = ' +
+            str((mp.P2.D/2)**2/(mp.d_dm1_dm2*mp.lambda0)))
 
     return out
 
@@ -171,7 +172,6 @@ def set_optional_variables(mp):
     if not hasattr(mp.path, 'jac'):
         mp.path.jac = os.path.join(mp.path.falco, 'data', 'jac')
 
-
     # Parallel processing
     if not hasattr(mp, "flagParallel"):
         mp.flagParallel = False
@@ -201,7 +201,7 @@ def set_optional_variables(mp):
         mp.est.ItrStartKF = 2  # Which iteration to start the Kalman filter at
     if not hasattr(mp.ctrl, 'flagUseModel'):
         mp.ctrl.flagUseModel = False  # Whether to perform a model-based (vs empirical) grid search for the controller
-    
+
     # Algorithmic Differentiation EFC options
     if not hasattr(mp.ctrl, 'ad'):
         mp.ctrl.ad = falco.config.Object()