Refactor done but sim results are broken

depr/shane_september_2023.ipynb

@@ -464,7 +464,7 @@
    "source": [
     "### Daytime Testing\n",
     "\n",
-    "How do we actually use this? This is where the `ShaneLantern` object comes in. The basic functions are `command_to_dm` and `get_image`. We keep track of the current DM command in `lantern.curr_dmc`, and can update it to put new patterns on."
+    "How do we actually use this? This is where the `ShaneLantern` object comes in. The basic functions are `command_to_dm` and `get_image`. We keep track of the current DM command in `lantern.actuators`, and can update it to put new patterns on."
    ]
   },
   {
@@ -572,7 +572,7 @@
       "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
       "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
       "Cell \u001b[0;32mIn[19], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m lantern\u001b[39m.\u001b[39;49mzern_to_dm(\u001b[39m2\u001b[39;49m, \u001b[39m1\u001b[39;49m)\n\u001b[1;32m      2\u001b[0m \u001b[39m# do I mean amp = 1.0 on Zernike 2, or 2.0 on Zernike 1?\u001b[39;00m\n",
-      "File \u001b[0;32m~/projects/photonics/photonics/shane.py:30\u001b[0m, in \u001b[0;36mShaneLantern.zern_to_dm\u001b[0;34m(self, z, amp)\u001b[0m\n\u001b[1;32m     28\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mzern_to_dm\u001b[39m(\u001b[39mself\u001b[39m, z, amp):\n\u001b[1;32m     29\u001b[0m     \u001b[39massert\u001b[39;00m \u001b[39mtype\u001b[39m(z) \u001b[39m==\u001b[39m \u001b[39mint\u001b[39m, \u001b[39m\"\u001b[39m\u001b[39mfirst argument must be an integer (Zernike number)\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m---> 30\u001b[0m     \u001b[39massert\u001b[39;00m \u001b[39mtype\u001b[39m(amp) \u001b[39m==\u001b[39m \u001b[39mfloat\u001b[39m, \u001b[39m\"\u001b[39m\u001b[39msecond argument must be a float (amplitude)\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m     31\u001b[0m     \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcurr_dmc[:] \u001b[39m=\u001b[39m \u001b[39m0.0\u001b[39m\n\u001b[1;32m     32\u001b[0m     \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcurr_dmc[z\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m] \u001b[39m=\u001b[39m amp\n",
+      "File \u001b[0;32m~/projects/photonics/photonics/shane.py:30\u001b[0m, in \u001b[0;36mShaneLantern.zern_to_dm\u001b[0;34m(self, z, amp)\u001b[0m\n\u001b[1;32m     28\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mzern_to_dm\u001b[39m(\u001b[39mself\u001b[39m, z, amp):\n\u001b[1;32m     29\u001b[0m     \u001b[39massert\u001b[39;00m \u001b[39mtype\u001b[39m(z) \u001b[39m==\u001b[39m \u001b[39mint\u001b[39m, \u001b[39m\"\u001b[39m\u001b[39mfirst argument must be an integer (Zernike number)\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m---> 30\u001b[0m     \u001b[39massert\u001b[39;00m \u001b[39mtype\u001b[39m(amp) \u001b[39m==\u001b[39m \u001b[39mfloat\u001b[39m, \u001b[39m\"\u001b[39m\u001b[39msecond argument must be a float (amplitude)\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m     31\u001b[0m     \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mactuators[:] \u001b[39m=\u001b[39m \u001b[39m0.0\u001b[39m\n\u001b[1;32m     32\u001b[0m     \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mactuators[z\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m] \u001b[39m=\u001b[39m amp\n",
       "\u001b[0;31mAssertionError\u001b[0m: second argument must be a float (amplitude)"
      ]
     }

photonics/experiments/deformable_mirrors.py

@@ -1,12 +1,14 @@
 import time
 import zmq
+from abc import ABC
 from tqdm import trange
 import numpy as np
 import hcipy as hc
 
 class ShaneDM:
     def __init__(self):
-        self.curr_dmc = np.zeros(Nmodes)
+        self.Nmodes = 12
+        self.actuators = np.zeros(self.Nmodes)
         self.setup_paramiko()
 
     def setup_paramiko(self):
@@ -24,28 +26,29 @@ def __del__(self):
     def apply_mode(self, z, amp, verbose=True):
         assert isinstance(z, int), "first argument must be an integer (Zernike number)"
         assert isinstance(amp, float), "second argument must be a float (amplitude)"
-        self.curr_dmc[:] = 0.0
-        self.curr_dmc[z-1] = amp
+        self.actuators[:] = 0.0
+        self.actuators[z-1] = amp
         self.command_to_dm(verbose=verbose)
 
-    def command_to_dm(self, verbose=True):
+    def command_to_dm(self, p, verbose=True):
         """
         Send a command to the ShaneAO woofer.
 
         Parameters:
             amplitudes - list or np.ndarray
             The amplitude of each mode in [1, 2, ..., Nmodes], in order.
         """
-        assert len(self.curr_dmc) == self.Nmodes, "wrong number of modes specified"
-        assert np.all(np.abs(self.curr_dmc) <= 1.0), "sending out-of-bounds amplitudes"
-        command = ",".join(map(str, self.curr_dmc))
+        self.actuators[:] = p[:]
+        assert len(self.actuators) == self.Nmodes, "wrong number of modes specified"
+        assert np.all(np.abs(self.actuators) <= 1.0), "sending out-of-bounds amplitudes"
+        command = ",".join(map(str, self.actuators))
         if verbose:
             print(f"DMC {command}.")
         #warnings.warn("If you see this and you're at Lick, uncomment the lines defining and running shell_command.")
         self.client.exec_command(f"/home/user/ShaneAO/shade/imageSharpen_nogui -s {command}")
 
     def send_zeros(self, verbose=True):
-        self.curr_dmc[:] = 0.0
+        self.actuators[:] = 0.0
         if verbose:
             print("Sending zeros.")
         self.command_to_dm(verbose=False)
@@ -90,6 +93,7 @@ def __del__(self):
 
 class SimulatedDM:
     def __init__(self, pupil_grid, dm_basis="modal", num_actuators=9, telescope_diameter=10):
+        self.Nmodes = num_actuators ** 2
         if dm_basis == "zonal":
             actuator_spacing = telescope_diameter / num_actuators
             influence_functions = hc.make_gaussian_influence_functions(pupil_grid, num_actuators, actuator_spacing)
@@ -105,6 +109,15 @@ def send_zeros(self, verbose=True):
             print("Sending zeros.")
         self.deformable_mirror.actuators[:] = 0.0
 
+    @property
+    def actuators(self):
+        # just for convenience; actuators setter is command_to_dm
+        # wait, actually, why don't I just...do that
+        return self.deformable_mirror.actuators
+
+    def command_to_dm(self, p, verbose=False):
+        self.deformable_mirror.actuators[:] = p[:]
+
     def apply_mode(self, mode, amplitude):
         self.deformable_mirror.actuators[mode] = amplitude
 

photonics/experiments/experiments.py

@@ -55,54 +55,49 @@ def measure_pl_flat(self):
 
     def openloop_experiment(self, patterns, tag="openloop_experiment", delay=0):
         start_stamp = datetime_ms_now()
-        self.send_zeros(verbose=False)
-        img = self.get_image()
+        self.dm.send_zeros(verbose=False)
+        img = self.wfs.get_image()
         start_time_stamp = time_ms_now()
-        bbox_shape = self.get_image().shape
+        bbox_shape = self.wfs.get_image().shape
         pl_readout = np.zeros((len(patterns), *bbox_shape))
 
         for (i, p) in enumerate(tqdm(patterns)):
-            self.curr_dmc = p
-            self.command_to_dm(verbose=False)
+            self.dm.command_to_dm(p, verbose=False)
             sleep(delay)
-            img = self.get_image()
+            img = self.wfs.get_image()
             pl_readout[i] = img
 
-        pl_intensities = np.array(list(map(self.get_intensities, pl_readout)))
-        self.send_zeros(verbose=False)
+        pl_intensities = np.array(list(map(self.wfs.get_intensities, pl_readout)))
+        self.dm.send_zeros(verbose=False)
         with h5py.File(self.filepath(f"{tag}_{datetime_now()}", ext="hdf5"), "w") as f:
             dmcs_dataset = f.create_dataset("dmc", data=patterns)
             if self.save_full_frame:
                 f.create_dataset("pl_images", data=pl_readout)
             pl_intensities_dataset = f.create_dataset("pl_intensities", data=pl_intensities)
-            pl_intensities_dataset.attrs["exp_ms"] = self.exp_ms
-            pl_intensities_dataset.attrs["gain"] = self.gain
-            pl_intensities_dataset.attrs["nframes"] = self.nframes
-            pl_intensities_dataset.attrs["centroids_dt"] = self.centroids_dt
+            pl_intensities_dataset.attrs["exp_ms"] = self.wfs.exp_ms
+            pl_intensities_dataset.attrs["gain"] = self.wfs.gain
+            pl_intensities_dataset.attrs["nframes"] = self.wfs.nframes
+            pl_intensities_dataset.attrs["centroids_dt"] = self.wfs.centroids_dt
 
-        self.destroy_paramiko()
         return pl_intensities
 
     def save_current_image(self, tag=""):
-        self.save(f"pl_image_{datetime_now()}_{tag}", self.get_image())
+        self.save(f"pl_image_{datetime_now()}_{tag}", self.wfs.get_image())
 
-    def make_interaction_matrix(self, amp_calib=0.01, thres=1/30, nm=None):
-        self.setup_paramiko()
-        if nm is None:
-            nm = self.Nmodes
-        self.int_mat = np.zeros((self.Nports, nm))
+    def make_interaction_matrix(self, nm, amp_calib=0.01, thres=1/30):
+        self.int_mat = np.zeros((self.wfs.Nports, nm))
         pushes = []
         pulls = []
         for i in trange(nm):
-            self.zern_to_dm(i + 1, amp_calib, verbose=False)
+            self.dm.apply_mode(i + 1, amp_calib)
             sleep(0.1)
-            s_push = self.get_image()
+            s_push = self.wfs.get_image()
             pushes.append(s_push)
-            self.zern_to_dm(i + 1, -amp_calib, verbose=False)
+            self.dm.apply_mode(i + 1, -amp_calib)
             sleep(0.1)
-            s_pull = self.get_image()
+            s_pull = self.wfs.get_image()
             pulls.append(s_pull)
-            s = (self.get_intensities(s_push) - self.get_intensities(s_pull)) / (2 * amp_calib)
+            s = (self.wfs.get_intensities(s_push) - self.wfs.get_intensities(s_pull)) / (2 * amp_calib)
             self.int_mat[:,i] = s.ravel()
 
         self.cmd_mat = np.linalg.pinv(self.int_mat, thres)
@@ -112,14 +107,13 @@ def make_interaction_matrix(self, amp_calib=0.01, thres=1/30, nm=None):
             pulls_dset = f.create_dataset("pull_frames", data=np.array(pulls))
             intmat_dset = f.create_dataset("intmat", data=self.int_mat)
             intmat_dset.attrs["amp_calib"] = amp_calib
-            intmat_dset.attrs["exp_ms"] = self.exp_ms
-            intmat_dset.attrs["gain"] = self.gain
-            intmat_dset.attrs["nframes"] = self.nframes
+            intmat_dset.attrs["exp_ms"] = self.wfs.exp_ms
+            intmat_dset.attrs["gain"] = self.wfs.gain
+            intmat_dset.attrs["nframes"] = self.wfs.nframes
             cmdmat_dset = f.create_dataset("cmdmat", data=self.cmd_mat)
             cmdmat_dset.attrs["thres"] = thres
 
-        self.send_zeros()
-        self.destroy_paramiko()
+        self.dm.send_zeros()
 
     def load_interaction_matrix(self, fname):
         with h5py.File(fname) as f:
@@ -129,43 +123,42 @@ def load_interaction_matrix(self, fname):
             self.cmd_mat = np.array(f["cmdmat"])
 
     def pseudo_cl_iteration(self, gain=0.1, verbose=True):
-        dm_start = np.copy(self.curr_dmc)
+        dm_start = np.copy(self.dm.actuators)
         if verbose:
-            print(f"Initial error = {rms(self.curr_dmc)}")
-        lantern_reading = np.zeros(self.Nmodes)
-        lantern_image = self.get_image()
-        lr = self.cmd_mat @ (self.get_intensities(lantern_image) - self.pl_flat)
+            print(f"Initial error = {rms(self.dm.actuators)}")
+        lantern_reading = np.zeros(self.dm.Nmodes)
+        lantern_image = self.wfs.get_image()
+        lr = self.cmd_mat @ (self.wfs.get_intensities(lantern_image) - self.pl_flat)
         lantern_reading[:len(lr)] = lr
         # does the sign of measured WF errors match the actual signs?
-        sign_match = ''.join(map(lambda x: str(int(x)), (np.sign(lantern_reading * self.curr_dmc) * 2 + 2)/4))
+        sign_match = ''.join(map(lambda x: str(int(x)), (np.sign(lantern_reading * self.dm.actuators) * 2 + 2)/4))
         if verbose:
             print(f"{sign_match}")
-        self.curr_dmc -= gain * lantern_reading
-        self.command_to_dm(verbose=verbose)
-        dm_end = np.copy(self.curr_dmc)
+        self.dm.command_to_dm(self.dm.actuators - gain * lantern_reading, verbose=verbose)
+        dm_end = np.copy(self.dm.actuators)
         if verbose:
-            print(f"Final error = {rms(self.curr_dmc)}")
+            print(f"Final error = {rms(self.dm.actuators)}")
 
         if not verbose:
             return lantern_image, lantern_reading
 
     def closed_loop(self, gain=0.1, niter=20):
-        dmcs = [np.copy(self.curr_dmc)]
+        dmcs = [np.copy(self.dm.actuators)]
         lantern_images, lantern_readings = [], []
         for i in trange(niter):
            lantern_image, lantern_reading = self.pseudo_cl_iteration(gain=gain, verbose=False)
-           dmcs.append(np.copy(self.curr_dmc))
+           dmcs.append(np.copy(self.dm.actuators))
            lantern_images.append(lantern_image)
            lantern_readings.append(lantern_reading)
 
         with h5py.File(self.filepath(f"closedloop_{datetime_now()}", ext="hdf5"), "w") as f:
             dmcs_dset = f.create_dataset("dmcs", data=np.array(dmcs))
-            dmcs_dset.attrs["exp_ms"] = self.exp_ms
-            dmcs_dset.attrs["gain"] = self.gain
-            dmcs_dset.attrs["centroids_timestamp"] = self.centroids_dt
+            dmcs_dset.attrs["exp_ms"] = self.wfs.exp_ms
+            dmcs_dset.attrs["gain"] = self.wfs.gain
+            dmcs_dset.attrs["centroids_timestamp"] = self.wfs.centroids_dt
             dmcs_dset.attrs["cmd_timestamp"] = self.cmd_dt
             images_dset = f.create_dataset("lantern_images", data=np.array(lantern_images))
-            images_dset.attrs["spot_radius_px"] = self.spot_radius_px
+            images_dset.attrs["spot_radius_px"] = self.wfs.spot_radius_px
             f.create_dataset("lantern_readings", data=np.array(lantern_readings))
 
         return dmcs, lantern_readings
@@ -174,14 +167,11 @@ def closed_loop_replay(self, fname):
         with h5py.File(fname) as f:
             dmcs = np.array(f["dmcs"])
         for (i, dmc) in enumerate(dmcs):
-            self.curr_dmc = dmc
-            self.command_to_dm()
+            self.dm.command_to_dm(dmc)
             # I don't have programmatic access to the PSF camera, so I'll just wait for user input
             input(f"Frame {i}.")
 
-    def measure_linearity(self, min_amp=-0.06, max_amp=0.06, step=0.01, nmodes=None):
-        if nmodes is None:
-            nmodes = self.Nmodes
+    def measure_linearity(self, nmodes, min_amp=-0.06, max_amp=0.06, step=0.01):
         all_recon = []
         amps = None
         for z in range(1, nmodes+1):
@@ -199,8 +189,8 @@ def measure_linearity(self, min_amp=-0.06, max_amp=0.06, step=0.01, nmodes=None)
 
     # different types of experiment
     def sweep_mode(self, z, min_amp=-1.0, max_amp=1.0, step=0.1):
-        amps = np.arange(min_amp, max_amp+2*step, step)
-        patterns = np.zeros((len(amps), self.Nmodes))
+        amps = np.arange(min_amp, max_amp+step/2, step)
+        patterns = np.zeros((len(amps), self.dm.Nmodes))
         patterns[:,z-1] = amps
         return amps, self.openloop_experiment(patterns, tag=f"sweep_mode_{z}")
 
@@ -212,7 +202,7 @@ def sweep_all_modes(self, min_amp=-1.0, max_amp=1.0, step=0.1):
         return amps, self.openloop_experiment(patterns, tag="sweep_all_modes")
 
     def random_combinations(self, Niters, lim=0.05):
-        inputzs = np.random.uniform(-lim, lim, (Niters+1, self.Nmodes))
+        inputzs = np.random.uniform(-lim, lim, (Niters+1, self.dm.Nmodes))
         self.openloop_experiment(inputzs, tag="random_combinations")
 
     def exptime_sweep(self, exp_ms_values, tag=""):
@@ -244,7 +234,7 @@ def snr_per_port(self, exp_ms_low, exp_ms_high, exp_ms_step):
         return exp_ms_values, snr
 
     def stability(self, n_pictures=20, wait_s=10):
-        patterns = np.tile(self.curr_dmc, (n_pictures,1))
+        patterns = np.tile(self.dm.actuators, (n_pictures,1))
         self.openloop_experiment(patterns, tag="stability", delay=wait_s)
 
     def stability_better(self, n_pictures=20, wait_s=2):

photonics/experiments/lantern_cameras.py

@@ -28,10 +28,9 @@ def measure_pl_flat(self):
         self.save(f"pl_flat_{datetime_now()}", self.pl_flat)
 
     def plot_ports(self, save=False):
-        # only run after set_centroids
         sc = plt.scatter(self.xc, self.yc, c=self.radial_shell)
-        plt.xlim((0, self.imgshape[1]))
-        plt.ylim((0, self.imgshape[0]))
+        plt.xlim((0, self.xg.shape[1]))
+        plt.ylim((0, self.xg.shape[0]))
         plt.xticks([])
         plt.yticks([])
         sc.axes.invert_yaxis()
@@ -50,10 +49,11 @@ def set_centroids(self, image, rerun=False):
             self.tag = input("enter a tag: ")
         centroids_file = os.path.join(DATA_PATH, "current_centroids", f"current_centroids_{self.tag}.hdf5")
         if os.path.exists(centroids_file) and not rerun:
-            with h5py.File(centroids_file, "w") as f:
+            with h5py.File(centroids_file, "r") as f:
                 centroids = np.array(f["centroids"])
                 self.xc = centroids[:,0]
                 self.yc = centroids[:,1]
+                self.radial_shell = centroids[:,2]
                 self.spot_radius_px = f["centroids"].attrs["spot_radius_px"]
             return
         global coords
@@ -76,7 +76,7 @@ def set_centroids(self, image, rerun=False):
 
             accepted_positions = input("Identified PL ports OK? [y/n] ") == "y"
 
-        self.xc, self.yc = ports_in_radial_order(self.refine_centroids(centroids, image))
+        self.xc, self.yc, self.radial_shell = ports_in_radial_order(self.refine_centroids(centroids, image))
         self.yi, self.xi = np.indices(image.shape)
         self.masks = [(self.xi - xc) ** 2 + (self.yi - yc) ** 2 <= self.spot_radius_px ** 2 for (xc, yc) in zip(self.xc, self.yc)]
         self.centroids_dt = datetime_now()
@@ -96,9 +96,10 @@ def set_centroids(self, image, rerun=False):
                 self.masks = [(self.xi - xc) ** 2 + (self.yi - yc) ** 2 <= self.spot_radius_px ** 2 for (xc, yc) in zip(self.xc, self.yc)]
             else:
                 accepted_radius = True
-        new_centroids = np.zeros_like(centroids)
+        new_centroids = np.zeros((len(centroids), 3))
         new_centroids[:,0] = self.xc
         new_centroids[:,1] = self.yc
+        new_centroids[:,2] = self.radial_shell
 
         with h5py.File(centroids_file, "w") as f:
             c = f.create_dataset("centroids", data=new_centroids)
@@ -192,6 +193,7 @@ def __init__(self, dm, optics, tag=""):
         self.optics = optics
         self.lantern_optics = LanternOptics(optics)
         self.exp_ms = 1e5
+        self.gain = 1
         self.nframes = 10
         self.detector = NoisyDetector(detector_grid=self.lantern_optics.focal_grid)
         self.dm = dm

photonics/utils.py

@@ -97,7 +97,7 @@ def ports_in_radial_order(points):
         prev_idx += nhull
         points = np.delete(points, v, axis=0)
 
-    return np.flip(xnew), np.flip(ynew)
+    return np.flip(xnew), np.flip(ynew), np.flip(radial_shell)
 
 def normalize(x):
     return x / np.sum(x)

scripts_py/experiment_sim_testing.py

@@ -1,12 +1,18 @@
-from photonics import SimulatedDM, SimulatedLanternCamera, Experiments
-from photonics.simulations.optics import Optics
 from hcipy import imshow_field
 from matplotlib import pyplot as plt
 
+from photonics import SimulatedDM, SimulatedLanternCamera, Experiments
+from photonics.simulations.optics import Optics
+from photonics.linearity import plot_linearity
+
 optics = Optics(lantern_fnumber=10) # so I flood the input/can see all the ports
 dm = SimulatedDM(optics.pupil_grid)
 pl = SimulatedLanternCamera(dm, optics, "spie24_sim_pl")
 exp = Experiments(dm, pl)
 
 pl.measure_dark()
-pl.set_centroids(pl.get_image())
+img = pl.get_image()
+pl.set_centroids(img, rerun=False) # set to True if we want to redo centroiding
+exp.measure_pl_flat()
+exp.make_interaction_matrix(9)
+amps, all_recon = exp.measure_linearity(3)