Merge branch 'develop'

Author: nrodlin

.github/workflows/gh-pages.yml

@@ -0,0 +1,49 @@
+name: Deploy GitHub Pages
+
+on:
+  push:
+    branches:
+      - master
+
+permissions:
+  contents: write
+
+jobs:
+  deploy:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Configure Git Credentials
+        run: |
+          git config user.name github-actions[bot]
+          git config user.email 41898282+github-actions[bot]@users.noreply.github.com
+
+      - uses: actions/setup-python@v5
+        with:
+          # Use standard (non-free-threaded) Python — cffi is incompatible with 3.13t
+          python-version: '3.12'
+
+      - name: Cache pip dependencies
+        uses: actions/cache@v4
+        with:
+          key: mkdocs-${{ hashFiles('requirements.txt') }}
+          path: ~/.cache/pip
+          restore-keys: |
+            mkdocs-
+
+      - name: Install documentation dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install mkdocs mkdocs-material mkdocs-jupyter mkdocstrings[python] ruff
+
+      - name: Install SAOS and its dependencies
+        run: |
+          # Install numpy, scipy, torch (CPU-only, lightweight) so mkdocstrings can import SAOS
+          pip install numpy scipy opencv-contrib-python joblib
+          pip install torch --index-url https://download.pytorch.org/whl/cpu
+          pip install -e . --no-deps
+
+      - name: Deploy to GitHub Pages
+        run: mkdocs gh-deploy --force
+

README.md

@@ -4,13 +4,31 @@
 
 SAOS modifies the internal architecture of OOPAO, which was heavily managed through the telescope class, to offer a framework centered on lines of sight, implemented via the **LightPath** class. This architectural shift enables easier and more powerful parallelization, essential for solar AO simulations where the Field of View (FoV) is divided into multiple sub-directions to physically introduce extended field effects into the simulation.
 
-SAOS is currently in the testing phase of the first release. All development is ongoing in the **tech-development** branch. We recommend waiting a little longer before using SAOS in your work!
+A key feature of SAOS is that it is a **Python-only** repository. All parallelization is managed via **joblib** and threads, leveraging the new **free-threaded Python 3.13** to enable extensive and straightforward parallelization. Pytorch is extensively used in the code for costly operations, being ready for an upgrade to a GPU version.
 
-A key feature of SAOS is that it is a **Python-only** repository. All parallelization is managed via **joblib** and threads, leveraging the new **free-threaded Python 3.13** to enable extensive and straightforward parallelization. Pytorch is extensively used in the code for costly operations, being ready to an upgrade to a GPU version.
+---
 
-# Guidelines
+## 📖 Documentation
 
-- If you are primarily a night-time AO user and do not require multi-conjugate AO (MCAO) or you need a pyramid WFS simulations, we recommend using **OOPAO** instead. See: [OOPAO AO4ELT7 Proceedings (Héritier et al. 2023)](https://hal.science/AO4ELT7/hal-04402878v1).
+Full documentation — installation guides, pedagogical tutorials, and API reference — is available at:
+
+**➡️ [https://nrodlin.github.io/SAOS/](https://nrodlin.github.io/SAOS/)**
+
+To build the documentation locally:
+```bash
+pip install mkdocs mkdocs-material mkdocs-jupyter mkdocstrings[python]
+mkdocs serve
+```
+
+## ⚙️ Installation
+
+For optimal performance, we strongly recommend using a free-threaded Python environment (e.g., Python 3.13.2t). 
+A quick installation snippet is provided below. For detailed instructions, please check the [Installation Guide](docs/installation.md).
+
+```bash
+git clone https://github.com/nrodlin/SAOS.git
+cd SAOS
+git checkout develop
+pip install -e .
+```
 
-- If you need to simulate multiple lines of sight or perform Solar AO simulations, **SAOS** is the repository for you.  
-  We have not yet published a full article describing SAOS, but it is in preparation —please keep an eye on the repository for updates!

SAOS/Atmosphere.py

@@ -19,7 +19,7 @@
 from joblib import Parallel, delayed
 
 import numpy as np
-import math
+
 import torch
 
 # Self dependencies
@@ -135,13 +135,11 @@ def __init__(self,
 
     def initializeAtmosphere(self, randomState=None):
         """
-        Initialize the atmosphere layers and associate them with a telescope.
+        Initialize the atmosphere layers using the telescope linked at construction time.
 
         Parameters
         ----------
-        telescope : Telescope or None
-            Telescope object to derive spatial and temporal resolution.
-        randomState : int or None
+        randomState : int or None, optional
             Seed for reproducible random number generation, by default None.
 
         Returns

SAOS/Controller.py

@@ -1,7 +1,7 @@
 import numpy as np
 import torch
 
-import h5py
+
 import time
 
 import logging
@@ -66,9 +66,12 @@ def __init__(self,
         # Default will change to list of size nDMs once the IM is scanned
         self.rcond = kwargs.get('rcond', 0.025)
         self.beta = kwargs.get('beta', 1e-4) # adim, adjusted through trial-error
+        
+        # Mask provided by the user to select specific WFS-DM links
+        self.control_mask = kwargs.get('control_mask', None)
 
         # Run the initialization of the reconstructor
-        self.reconstructor, self.modal_basis, self.mask, self.altitude = self.initializeReconstructor(self.reconstructionMethod, interactionMatrix)                
+        self.reconstructor, self.modal_basis, self.mask, self.discarded_modes = self.initializeReconstructor(self.reconstructionMethod, interactionMatrix)                
 
         # Setup the controller
 
@@ -103,6 +106,27 @@ def __init__(self,
         self.initializeController(self.controllerType, self.reconstructor)
 
     def initializeReconstructor(self, reconstructionMethod, interactionMatrix):
+        """
+        Initialize the reconstructor matrix from the measured interaction matrices.
+
+        Parameters
+        ----------
+        reconstructionMethod : str
+            Type of reconstructor ('inversion' or 'tikhonov').
+        interactionMatrix : InteractionMatrixHandler
+            Object containing the measured interaction matrices and modal basis.
+
+        Returns
+        -------
+        reconstructor : list
+            List of reconstructor matrices per DM.
+        modal_basis : list
+            List of modal basis per DM.
+        mask : np.ndarray
+            Boolean mask indicating interactions between DMs and light paths.
+        discarded_modes : list
+            List of number of discarded modes per DM.
+        """
         self.logger.info('Controller::initializeReconstructor - Computing the reconstructor.')
         t0 = time.time()
 
@@ -127,6 +151,20 @@ def initializeReconstructor(self, reconstructionMethod, interactionMatrix):
                 if interactionMatrix.interaction_matrix_warehouse[i][j]['IM'] is not None:
                     mask[i, j] = True
 
+        if hasattr(self, 'control_mask') and self.control_mask is not None:
+            # Check dimensions
+            control_mask_arr = np.array(self.control_mask, dtype=bool)
+            if control_mask_arr.shape != (nDMs, nLPs):
+                self.logger.error(f'Controller - control_mask shape must be ({nDMs}, {nLPs})')
+                raise ValueError(f'control_mask shape mismatch. Expected ({nDMs}, {nLPs}), got {control_mask_arr.shape}')
+            
+            # Warn if user requests control where no IM exists
+            invalid_requests = control_mask_arr & (~mask)
+            if np.any(invalid_requests):
+                self.logger.warning('Controller - control_mask requests control for DM/LP pairs without an interaction matrix. These will be ignored.')
+                
+            mask = mask & control_mask_arr
+
         # Check the reconstructor parameters
         if reconstructionMethod == 'inversion':
             if isinstance(self.rcond, list):
@@ -154,10 +192,15 @@ def initializeReconstructor(self, reconstructionMethod, interactionMatrix):
                     modal_basis_type = interactionMatrix.interaction_matrix_warehouse[i][j]['modalBasis']
                     modal_basis.append(torch.as_tensor(interactionMatrix.modal_basis[i][modal_basis_type], dtype=torch.float64, device=self.device))
                     break
-        # Get altitudes:
-        altitude = []
-        for i in range(len(interactionMatrix.dm_scanned_list)):
-            altitude.append(interactionMatrix.dm_scanned_list[i].altitude)
+        # Get discarded modes metadata:
+        discarded_modes = []
+        for i in range(nDMs):
+            found_discarded_modes = 0
+            for j in range(nLPs):
+                if interactionMatrix.interaction_matrix_warehouse[i][j]['IM'] is not None:
+                    found_discarded_modes = interactionMatrix.interaction_matrix_warehouse[i][j].get('discarded_modes', 0)
+                    break
+            discarded_modes.append(found_discarded_modes)
 
         # Now, define the reconstruction matrices for each DM
 
@@ -170,26 +213,48 @@ def initializeReconstructor(self, reconstructionMethod, interactionMatrix):
                     # Append the IMs to shape one large matrix of size nValidAct x nSignals
                     interaction_matrix_per_DM.append(interactionMatrix.interaction_matrix_warehouse[i][j]['IM'])
             # Compute the reconstructor
-            interaction_matrix_per_DM = torch.as_tensor(np.vstack(interaction_matrix_per_DM), dtype=torch.float64, device=self.device).squeeze()
-            if reconstructionMethod == 'inversion':
-                temp_reconstructor = torch.linalg.pinv(interaction_matrix_per_DM, self.rcond[i])
-            elif reconstructionMethod == 'tikhonov':
-                # (D.T@D + alfa*I)@D.T --> implemented through SVD to improve the stability of the inversion and the automation of alfa
-                H = interaction_matrix_per_DM
-                U, S, Vh = torch.linalg.svd(H, full_matrices=False)
-                alfa = self.beta[i] * torch.max(S)**2
-                S_reg = S / (S**2 + alfa)
-                temp_reconstructor = Vh.T @ torch.diag(S_reg) @ U.T
+            if len(interaction_matrix_per_DM) == 0:
+                self.logger.warning(f'Controller - DM {i} has no associated WFS in the control mask. Setting reconstructor to zero.')
+                nModes = modal_basis[i].shape[1]
+                temp_reconstructor = torch.zeros((nModes, 0), dtype=torch.float64, device=self.device)
             else:
-                self.logger.error('Controller::initializeReconstructor - Unknown reconstructor')
-                raise ValueError('Unknown reconstructor method.')
+                interaction_matrix_per_DM = torch.as_tensor(np.vstack(interaction_matrix_per_DM), dtype=torch.float64, device=self.device).squeeze()
+                if interaction_matrix_per_DM.ndim == 1:
+                    interaction_matrix_per_DM = interaction_matrix_per_DM.unsqueeze(0)
+                if reconstructionMethod == 'inversion':
+                    temp_reconstructor = torch.linalg.pinv(interaction_matrix_per_DM, self.rcond[i])
+                elif reconstructionMethod == 'tikhonov':
+                    # (D.T@D + alfa*I)@D.T --> implemented through SVD to improve the stability of the inversion and the automation of alfa
+                    H = interaction_matrix_per_DM
+                    U, S, Vh = torch.linalg.svd(H, full_matrices=False)
+                    alfa = self.beta[i] * torch.max(S)**2
+                    S_reg = S / (S**2 + alfa)
+                    temp_reconstructor = Vh.T @ torch.diag(S_reg) @ U.T
+                else:
+                    self.logger.error('Controller::initializeReconstructor - Unknown reconstructor')
+                    raise ValueError('Unknown reconstructor method.')
             reconstructor.append(temp_reconstructor)
 
         self.logger.info(f'Controller::initializeReconstructor - Reconstruction took {time.time()-t0}[s]')
 
-        return reconstructor, modal_basis, mask, altitude
+        return reconstructor, modal_basis, mask, discarded_modes
 
     def initializeController(self, controllerType, reconstructor):
+        """
+        Initialize the control state (history buffers) based on the controller type.
+
+        Parameters
+        ----------
+        controllerType : str
+            Type of controller ('leaky', 'forwardPI', 'backwardPI').
+        reconstructor : list
+            List of reconstructor matrices per DM.
+
+        Returns
+        -------
+        bool
+            True if initialization succeeds.
+        """
 
         if controllerType == 'leaky':
             self.command_previous = [torch.zeros((reconstructor[i].shape[0],1), dtype=torch.float64, device=self.device) for i in range(len(reconstructor))]
@@ -202,6 +267,19 @@ def initializeController(self, controllerType, reconstructor):
         return True
 
     def computeControlAction(self, lightPaths):
+        """
+        Compute the control action for each DM given the wavefront error from the light paths.
+
+        Parameters
+        ----------
+        lightPaths : list
+            List of LightPath objects that contain the wavefront error measurements.
+
+        Returns
+        -------
+        dm_cmd : list
+            List of command arrays to be sent to each Deformable Mirror.
+        """
 
         # Get the combined measurement array for each DM
         error = []
@@ -213,7 +291,10 @@ def computeControlAction(self, lightPaths):
                     combined_slopes.append(lightPaths[j].get_wavefront_error())
 
             # Convert to torch
-            error.append((-1)*torch.as_tensor(np.hstack(combined_slopes).T, dtype=torch.float64, device=self.device).unsqueeze(1)) # -1 for the feedback
+            if len(combined_slopes) > 0:
+                error.append((-1)*torch.as_tensor(np.hstack(combined_slopes).T, dtype=torch.float64, device=self.device).unsqueeze(1)) # -1 for the feedback
+            else:
+                error.append(torch.zeros((0, 1), dtype=torch.float64, device=self.device))
 
         # Compute the DM command
         modal_error = []
@@ -235,10 +316,8 @@ def computeControlAction(self, lightPaths):
         dm_cmd = []
 
         for i in range(len(self.reconstructor)):
-            if self.altitude[i] > 0: # TT is discarded automatically in the IM measurement
-                dm_cmd.append(self.modal_basis[i][:,2:2+self.reconstructor[i].shape[0]] @ modal_cmd[i])
-            else:
-                dm_cmd.append(self.modal_basis[i][:,:self.reconstructor[i].shape[0]] @ modal_cmd[i])
+            offset = self.discarded_modes[i]
+            dm_cmd.append(self.modal_basis[i][:, offset : offset + self.reconstructor[i].shape[0]] @ modal_cmd[i])
 
         # Update history buffers for the next iteration
 

SAOS/CorrelatingShackHartmann.py

@@ -12,7 +12,7 @@
 from queue import Queue
 
 import numpy as np
-import scipy as sp
+
 import torch
 
 import cv2
@@ -67,8 +67,6 @@ def __init__(self,
             Extra FoV in [arcsec] that is taken for the FFT computation, in order to reduce wrapping effects.
         use_brightest : int, optional
             Picks the n brightest pixels as threshold for center-of-gravity spot detection.
-        is_geometric : bool, optional
-            Enable geometric mode (gradient-based measurement).
         threshold_convolution : float, optional
             Cut-off threshold for Gaussian convolution.
         unit_in_rad : bool, optional
@@ -568,6 +566,23 @@ def get_psf(self, phase, fwhm, npix_sun):
         return psf
 
     def compute_images(self, psf, subDirs_sun, new_px):
+        """
+        Compute correlated subaperture images using convolution.
+
+        Parameters
+        ----------
+        psf : torch.Tensor
+            PSF of the subapertures.
+        subDirs_sun : np.ndarray
+            Sun subdirectories used as object.
+        new_px : int
+            New size in pixels for the sun subdirectories.
+
+        Returns
+        -------
+        torch.Tensor
+            Convoluted images.
+        """
 
         # Convert to Tensor
         sun_torch = torch.from_numpy(subDirs_sun).contiguous().to(self.device) # 4D: 
@@ -601,6 +616,21 @@ def compute_images(self, psf, subDirs_sun, new_px):
         return sun_patches
 
     def merges_images(self, sun_patches, src):
+        """
+        Merge sun patches into a single frame.
+
+        Parameters
+        ----------
+        sun_patches : torch.Tensor
+            The convoluted sun patches per subaperture.
+        src : Source
+            The solar source object.
+
+        Returns
+        -------
+        np.ndarray
+            The merged full-field images.
+        """
         # Resize the 2D filter
         filter_2D_torch = torch.from_numpy(src.filter_2D).contiguous().to(self.device)
         filter_2D_torch = filter_2D_torch.view(filter_2D_torch.shape[0], filter_2D_torch.shape[1], -1).permute(2, 0, 1)
@@ -844,12 +874,11 @@ def wfs_measure(self,phase_in, src, pseudoref=None, reference_slopes=None):
             Phase map input [radians].
         src : Source
             Source object.
-        integrate : bool, optional
-            Whether to include camera integration effects.
         pseudoref : torch.Tensor
             Pseudo-reference image used for the correlations, if None the method picks one.
         reference_slopes : np.ndarray
             Reference for the WF sensor, if None, geometric references are used.
+
         Returns
         -------
         tuple