<jupyter_start><jupyter_text>FPM SimulationThis ntoebook is not fully functional.<jupyter_code>%matplotlib widget
%load_ext autoreload
%autoreload 2

import numpy as np
import matplotlib.pyplot as plt
import llops as yp
import comptic
import llops.operators as ops
from llops.display import progress_bar

# Import FPM module
from comptic import fpm

# Set default datatype
yp.setDefaultDatatype('complex32')
yp.setDefaultBackend('numpy')<jupyter_output>The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload<jupyter_text>Generate Object<jupyter_code>object_shape = (256, 256)
object_true = yp.simulation.ucb(object_shape)

plt.figure()
plt.subplot(121)
plt.imshow(yp.abs(object_true))
plt.title('Amplitude')
plt.subplot(122)
plt.imshow(yp.angle(object_true))
plt.title('Phase')<jupyter_output><empty_output><jupyter_text>Define Settings<jupyter_code># Generate system parameters
system_parameters = fpm.getDefaultSystemParams(pixel_count=object_shape)

# Show led positions
fpm.plotLedPattern( **system_parameters)<jupyter_output><empty_output><jupyter_text>Generate Pupil Function<jupyter_code># Generate Pupil
pupil = comptic.imaging.pupil(object_shape, **system_parameters)<jupyter_output><empty_output><jupyter_text>Generate LED Sequences<jupyter_code># Select LED pattern. Should be a list of list of integers; each internal list has the LED id positions which are illuminated.
single = [[i] for i, x in enumerate(system_parameters['illumination_source_position_list_na'])]
multiplexed = [[4 * i + j for j in range(4)] for i in range(8)]
patterns = single[:120]<jupyter_output><empty_output><jupyter_text>Generate Forward Models<jupyter_code># Generate Intensity Operator
I = ops.Intensity(object_shape)

# Generate Phase Ramp Operator
R = ops.PhaseRamp(object_shape)

# Generate Fourier Transform Operator
F = ops.FourierTransform(object_shape)

# Generate pupil dot multiplication operator
P = ops.Diagonalize(pupil, label='\\tilde P')

# Parse Illumination wavelength
illumination_wavelength = system_parameters['illumination_wavelength']

# Convert led position list to correct backend
led_position_list_na = yp.asarray(system_parameters['illumination_source_position_list_na'], dtype='complex32')

# Derived constants
effective_pixel_size = system_parameters['camera_pixel_size'] / system_parameters['objective_magnification'] 
dkxy = yp.asarray([1 / (sz * effective_pixel_size) for sz in object_shape], dtype='complex32')

# Initialise forward model list (A_list) and measurement list (y_list)
A_list, y_list = [], []
for pattern_index, led_pattern in enumerate(tqdm(patterns, 'Forward Models Generated')):
    A_pattern = None
    
    # Loop over all LEDs in this pattern
    for led_id in led_pattern:
        led_position_na = yp.squeeze(led_position_list_na[led_id])
        
        # Determine the shift of the LED position in pixels
        pixel_shift = led_position_na / illumination_wavelength / dkxy
        
        # Create phase ramp multiplication operator
        S = ops.Diagonalize(R * pixel_shift, label='S_{' + str(pattern_index) + '}')
        
        # Generate complete forward model
        _A = I * (F.H * P * F * S)
        
        # Add this forward operator to the list
        A_pattern = _A if not A_pattern else A_pattern + _A
    
    # Append forward model to list
    A_list.append(A_pattern)
    
    # Generate measurement
    y_list.append(A_pattern * object_true)
    
# Combine these operators into a set if vertically stacked measurements
A = ops.Vstack(A_list)

# Create measurements
y = ops.VecStack(y_list)<jupyter_output>Forward Models Generated: 100%|██████████| 120/120 [00:01<00:00, 66.36it/s]<jupyter_text>Create Objective Function<jupyter_code>objective = yp.solvers.objectivefunctions.L2(A, y) + yp.solvers.regularizers.WaveletSparsity(object_shape)

# Render Latex
objective.latex()
# Define Initialization
initialization = yp.ones(object_shape)

# Solve by accelerated gradient descent
x_opt = yp.solvers.Fista(objective).solve(iteration_count=20, display_type='text', step_size=5e-3,
                                          use_nesterov_acceleration=True, nesterov_restart_enabled=True, initialization=initialization)

# Show recovered object
plt.figure()
plt.subplot(121)
plt.imshow(np.abs(x_opt), cmap = 'gray')
plt.title('Reconstructed Amplitude')

plt.subplot(122)
plt.imshow(np.angle(x_opt), cmap = 'gray')
plt.title('Reconstructed Phase')

plt.tight_layout()
plt.show()<jupyter_output>Minimizing function:
StarcoderdataJupyterStructuredCleanDedup