<jupyter_start><jupyter_text>Drifting coro floorIncorporating a drifting coro floor into the calculation of the PASTIS matrix  MotivationThroughout all our previous work, the coronagraph floor (contrast floor) $c_0$ has been assumed to be constant across the time we need to measure the full PASTIS matrix. As we have seen on HiCAT though, this is not the case, so we need to allow for a drifting contrast floor in the semi-analytic calculation of the PASTIS matrix. Reminder: static coro floorIn the theoretical JATIS 2021 paper, we have Eq. 15:$$c_{ij} = c_0 + a_c^2 m_{ii} + a_c^2 m_{jj} + 2 a_c^2 m_{ij}$$Solving for the diagonal elements $m_{ii}$ yields ($m_{ij} = 0$):$$m_{ii} = \frac{c_{ii} - c_0}{a_c^2}$$Using the above and solving for $m_{ij}$ yields:$$m_{ij} = \frac{c_{ij} + c_0 - c_{ii} - c_{jj}}{2 a_c^2}$$ Development for drifting coro floorThe above assumes that $c_0$ stays constant across all measurements of $c_{ij}$. This is not true on HiCAT, so we have to work in a $c_0$ that depends on $c_{ij}$. We can do this starting from the Eq. 15 in the JATIS paper (see above). The equation for the diagonal elements $m_{ii}$ stays the same, but we now want to derive the non-diagonal elements from the already calculated diagonal elements, not the $c_{ij}$. In that way, each measurement $c_{ij}$, and each matrix element $m_{ij}$ only depends on its specific, time-dependent measurement of $c_0$.In this case:$$m_{ii} = \frac{c_{ii} - c_{0_{ii}}}{a_c^2}$$and:$$m_{ij} = \frac{c_{ij} - c_{0_{ij}}}{2 a_c^2} - \frac{m_{ii}}{2} - \frac{m_{jj}}{2}$$  1: Making sure constant $c_0$ still works after refactor Comparing some matrix results<jupyter_code>import os
from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np

from pastis.matrix_building_numerical import pastis_from_contrast_matrix
import pastis.util
fname_old = '/Users/ilaginja/data_from_repos/pastis_data/2021-01-08T23-22-15_luvoir-small_develop/matrix_numerical/PASTISmatrix_num_piston_Noll1.fits'
fname_new = '/Users/ilaginja/data_from_repos/pastis_data/2021-01-09T01-01-53_luvoir-small_drift_0842c7/matrix_numerical/PASTISmatrix_num_piston_Noll1.fits'
pmatrix_old = fits.getdata(fname_old)
pmatrix_new = fits.getdata(fname_new)
plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
plt.imshow(pmatrix_old, origin='lower')
plt.title('pmatrix_old')
plt.colorbar()

plt.subplot(1, 3, 2)
plt.imshow(pmatrix_new, origin='lower')
plt.title('pmatrix_new')
plt.colorbar()

plt.subplot(1, 3, 3)
plt.imshow(pmatrix_old - pmatrix_new, origin='lower')
plt.title('diff')

plt.colorbar()
print(pmatrix_old[50,50])
print(pmatrix_new[50,50])
print(pmatrix_old[50,50] - pmatrix_new[50,50])
print(pmatrix_old[45,111])
print(pmatrix_new[45,111])
print(pmatrix_old[45,111] - pmatrix_new[45,111])
a = np.array([[1,2,3,],[4,5,6],[7,8,9]])
b = np.array([[1,24,3,],[4,15,6],[7,28,9]])
print(a)
print(b)
print(np.diag(b))
np.fill_diagonal(a, np.diag(b))
print(a)<jupyter_output><empty_output><jupyter_text>Compare contrast matrices<jupyter_code>cname_new = '/Users/ilaginja/data_from_repos/pastis_data/2021-01-08T13-02-27_luvoir-small/matrix_numerical/contrast_matrix.fits'
cname_old = '/Users/ilaginja/data_from_repos/pastis_data/2021-01-08T14-47-37_luvoir-small/matrix_numerical/pair-wise_contrasts.fits'
cname_new_norm_at_end = '/Users/ilaginja/data_from_repos/pastis_data/2021-01-08T19-28-09_luvoir-small/matrix_numerical/contrast_matrix.fits'
cm_new = fits.getdata(cname_new)
cm_old_subtracted = fits.getdata(cname_old)
cm_new_end = fits.getdata(cname_new_norm_at_end)
# cm_old had the coro floor subtracted, so I have to readd it - but only to the filled half of the matrix
coro_floor_old = 4.315823935036038e-11
coro_floor_matrix = np.zeros((120, 120))
seg_combos = list(pastis.util.segment_pairs_non_repeating(120))
len(seg_combos)
for pair in pastis.util.segment_pairs_non_repeating(120):    # this util function returns a generator
    if pair[0] != pair[1]:    # exclude diagonal elements
        coro_floor_matrix[pair[0], pair[1]] = coro_floor_old
# Fill diagonal
np.fill_diagonal(coro_floor_matrix, coro_floor_old)
cm_old = cm_old_subtracted + coro_floor_matrix
plt.figure(figsize=(10, 10))

plt.subplot(2, 2, 1)
plt.imshow(cm_new, origin='lower')
plt.title('cm_new')
plt.colorbar()

plt.subplot(2, 2, 2)
plt.imshow(cm_old, origin='lower')
plt.title('cm_old')
plt.colorbar()

plt.subplot(2, 2, 3)
plt.imshow(cm_new_end, origin='lower')
plt.title('cm_old_end')
plt.colorbar()

plt.subplot(2, 2, 4)
plt.imshow(cm_old-cm_new_end, origin='lower')
plt.title('diff')
plt.colorbar()
print(cm_old[50,50])
print(cm_new_end[50,50])
print(cm_old[50,50] - cm_new_end[50,50])
print(np.min(cm_old-cm_new_end))
# Top left triangle
print(cm_old[111,45])
print(cm_new_end[111,45])
print(cm_old[111,45] - cm_new_end[111,45])
# Bottom right triangle
print(cm_old[45,111])
print(cm_new_end[45,111])
print(cm_old[45,111] - cm_new_end[45,111])<jupyter_output><empty_output><jupyter_text>The contrat matrices are exactly the same, which is good. Now I can start debugging the analytical calculation procedure for the PASTIS matrix.I'll have a brief look at the difference anyway.<jupyter_code>diff = cm_old - cm_new
plt.figure(figsize=(10,10))
plt.imshow(diff)
plt.xlim(95,119)
plt.ylim(85, 119)
plt.colorbar()
print(diff[100,100])
print(diff[114,114])<jupyter_output><empty_output><jupyter_text>Calculate $M$ from $C$<jupyter_code>seglist = pastis.util.get_segment_list('LUVOIR')

# need uubtracted contrast matrices
#contrast_matrix = cm_old_subtracted
contrast_matrix = cm_new - coro_floor_matrix
#contrast_matrix = cm_new_end - coro_floor_matrix
## OLD WAY SHORT
# make sure you're at the old commit
matrix_pastis_dev_short = pastis_from_contrast_matrix(contrast_matrix, seglist, 1e-9)
plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
plt.imshow(matrix_pastis_dev_short, origin='lower')
plt.title('matrix_pastis_short')
plt.colorbar()

plt.subplot(1, 3, 2)
plt.imshow(pmatrix_old, origin='lower')
plt.title('pmatrix_old')
plt.colorbar()

plt.subplot(1, 3, 3)
plt.imshow(matrix_pastis_dev_short - pmatrix_old, origin='lower')
plt.title('diff')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>Difference is zero, as it should be. Good, moving on.<jupyter_code>## OLD WAY
# make sure you're at the old commit

# Create future (half filled) PASTIS matrix
matrix_pastis_half = np.copy(contrast_matrix)     # This will be the final PASTIS matrix.

# Calculate the off-axis elements in the (half) PASTIS matrix
for pair in pastis.util.segment_pairs_non_repeating(contrast_matrix.shape[0]):    # this util function returns a generator
    if pair[0] != pair[1]:    # exclude diagonal elements
        matrix_off_val = (contrast_matrix[pair[0], pair[1]] - contrast_matrix[pair[0], pair[0]] - contrast_matrix[pair[1], pair[1]]) / 2.
        matrix_pastis_half[pair[0], pair[1]] = matrix_off_val
        print(f'Off-axis for i{seglist[pair[0]]}-j{seglist[pair[1]]}: {matrix_off_val}')
# Symmetrize the half-PASTIS matrix
print('Symmetrizing PASTIS matrix')
matrix_pastis_dev = pastis.util.symmetrize(matrix_pastis_half)
print('Normalizing PASTIS matrix')
matrix_pastis_dev /= np.square(1e-9 * 1e9)
plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
plt.imshow(matrix_pastis_dev, origin='lower')
plt.colorbar()

plt.subplot(1, 3, 2)
plt.imshow(pmatrix_old, origin='lower')
plt.colorbar()

plt.subplot(1, 3, 3)
plt.imshow(matrix_pastis_dev - pmatrix_old, origin='lower')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>Ok this means I can reporoduce the old way step by step, both with the old as well as with the new contrast matrix. Whether the normalization happens at the beginning or at the end does not change the result.<jupyter_code>seglist = pastis.util.get_segment_list('LUVOIR')
wfe_aber = 1e-9

# need unsubtracted contrast matrices
contrast_matrix = cm_old
#contrast_matrix = cm_new #- coro_floor_matrix
#contrast_matrix = cm_new_end
# NEW WAY SHORT
# make sure you're at the new commit
matrix_pastis_commit = pastis_from_contrast_matrix(contrast_matrix, seglist, wfe_aber, coro_floor_old)
plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
plt.imshow(matrix_pastis_commit, origin='lower')
plt.colorbar()

plt.subplot(1, 3, 2)
plt.imshow(pmatrix_old, origin='lower')
plt.colorbar()

plt.subplot(1, 3, 3)
plt.imshow(matrix_pastis_commit - pmatrix_old, origin='lower')
plt.colorbar()
# NEW WAY
# Normalization
print('Normalization')
contrast_matrix /= np.square(wfe_aber * 1e9)  # 1e9 converts the calibration aberration back to nanometers
coro_floor = coro_floor_old / np.square(wfe_aber * 1e9)
# Create future (half filled) PASTIS matrix
matrix_pastis_half = np.zeros_like(contrast_matrix)     # This will be the final PASTIS matrix.
#matrix_pastis_half = np.copy(contrast_matrix)

# First calculate the on-axis elements, which just need to have the coronagraph floor subtracted
np.fill_diagonal(matrix_pastis_half, np.diag(contrast_matrix)-coro_floor)
#log.info('On-axis elements of PASTIS matrix calculated')

# Calculate the off-axis elements in the (half) PASTIS matrix
for pair in pastis.util.segment_pairs_non_repeating(contrast_matrix.shape[0]):    # this util function returns a generator
    if pair[0] != pair[1]:    # exclude diagonal elements
        matrix_off_val = (contrast_matrix[pair[0], pair[1]] + coro_floor - contrast_matrix[pair[0], pair[0]] - contrast_matrix[pair[1], pair[1]]) / 2.
        matrix_pastis_half[pair[0], pair[1]] = matrix_off_val
        #print(f'Off-axis for i{seglist[pair[0]]}-j{seglist[pair[1]]}: {matrix_off_val}')
# Symmetrize the half-PASTIS matrix
print('Symmetrizing PASTIS matrix')
matrix_pastis_commit = pastis.util.symmetrize(matrix_pastis_half)
# Old off-axis calculation

# Create future (half filled) PASTIS matrix
matrix_pastis_half = np.copy(contrast_matrix)    # This will be the final PASTIS matrix.

# Calculate the off-axis elements in the (half) PASTIS matrix
for pair in pastis.util.segment_pairs_non_repeating(contrast_matrix.shape[0]):    # this util function returns a generator
    if pair[0] != pair[1]:    # exclude diagonal elements
        matrix_off_val = (contrast_matrix[pair[0], pair[1]] - contrast_matrix[pair[0], pair[0]] - contrast_matrix[pair[1], pair[1]]) / 2.
        matrix_pastis_half[pair[0], pair[1]] = matrix_off_val
        print(f'Off-axis for i{seglist[pair[0]]}-j{seglist[pair[1]]}: {matrix_off_val}')
plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
plt.imshow(matrix_pastis_commit, origin='lower')
plt.colorbar()

plt.subplot(1, 3, 2)
plt.imshow(pmatrix_old, origin='lower')
plt.colorbar()

plt.subplot(1, 3, 3)
plt.imshow(matrix_pastis_commit - pmatrix_old, origin='lower')
plt.colorbar()
matrix_pastis_half = np.zeros_like(contrast_matrix)
np.fill_diagonal(matrix_pastis_half, np.diag(contrast_matrix))
plt.imshow(matrix_pastis_half)
plt.colorbar()
zer = np.zeros_like(contrast_matrix)
plt.imshow(zer)
plt.colorbar()
zer -= coro_floor_old
plt.imshow(zer)
plt.colorbar()
zer = np.triu(zer)
plt.imshow(zer)
plt.colorbar()<jupyter_output><empty_output><jupyter_text>2: Making sure it works for a drifting $c_0$To make sure this works well, I will use one of the contrast matrices above (they're all the same). I will generate a contrast floor matrix containing random values, remembering it needs to only be the "upper" triangle (`np.triu()`) and the rest zero. Then I will see whether the code produces a correct PASTIS matrix.<jupyter_code>coro_floor_matrix_full = np.random.normal(loc=0, scale=1, size=(120,120)) * coro_floor_old * 10
plt.imshow(coro_floor_matrix_full)
plt.colorbar()
coro_floor_matrix = np.triu(coro_floor_matrix_full)
plt.imshow(coro_floor_matrix)
plt.title('Random coro floor values per measurement')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>To create a contrast matrix that uses the above random coronagraph floor, I will do:- take a contrast matrix from which its original $c_0$ was already subtracted- add the random $c_0$ matrix to this subtracted contrast matrix- run it through the analytical PASTIS matrix generation that uses a drifting $c_0$- make sure it yields the same PASTIS matrix like this contrast matrix yielded when it had a constant $c_0$<jupyter_code># contrast matrix of choice:
contrast_matrix = cm_old_subtracted
plt.imshow(contrast_matrix)
plt.title('Original contrast matrix, no $c_0$')
plt.colorbar()
# add random c0 array
contrast_matrix += coro_floor_matrix
plt.imshow(contrast_matrix)
plt.title('Contrast matrix, with added random $c_0$')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>Calculate the PASTIS matrix from this.<jupyter_code>seglist = pastis.util.get_segment_list('LUVOIR')
wfe_aber = 1e-9
matrix_pastis_drift = pastis_from_contrast_matrix(contrast_matrix, seglist, wfe_aber, coro_floor_matrix)
diff = matrix_pastis_drift - pmatrix_old

plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
plt.imshow(matrix_pastis_drift, origin='lower')
plt.title('matrix_pastis_short')
plt.colorbar()

plt.subplot(1, 3, 2)
plt.imshow(pmatrix_old, origin='lower')
plt.title('pmatrix_old')
plt.colorbar()

plt.subplot(1, 3, 3)
plt.imshow(matrix_pastis_drift - pmatrix_old, origin='lower')
plt.title('diff')
plt.colorbar()
assert np.allclose(matrix_pastis_drift, pmatrix_old, rtol=1e-8, atol=1e-24), 'Calculated LUVOIR small PASTIS matrix is wrong.'
diff[40,40]<jupyter_output><empty_output><jupyter_text>Works.<jupyter_code># Make empty PASTIS matrix
matrix_pastis_half = np.zeros_like(contrast_matrix)
plt.imshow(matrix_pastis_half)
plt.colorbar()
# Fill diagonal elements
np.fill_diagonal(matrix_pastis_half, np.diag(contrast_matrix) - np.diag(coro_floor_matrix))
plt.imshow(matrix_pastis_half)
plt.colorbar()
assert (np.diag(matrix_pastis_half) == np.diag(pmatrix_old)).all, 'something wrong'<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12480940971465912, 'fraction_numerical': 0.033108255282073625, 'mean_word_length': 6.510359869138495, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '324115', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5517, 'n_tokens_neox': 5030, 'n_words': 1165, 'path': 'notebooks/LUVOIR-A/20_drifting coro floor.ipynb', 'repo_name': 'spacetelescope/PASTIS', 'size': '13773'}