desihub / specsim

Quick simulations of spectrograph response
2 stars 10 forks source link

Correct output resolution matrix #128

Closed p-slash closed 1 year ago

p-slash commented 1 year ago

This PR addresses issue #126, where the output resolution matrix from camera.py was averaged over rows and columns. This double downsampling gives incorrect resolution matrix outputs that was used in P1D mocks from quickspectra & quickquasars. The solution to eliminate one of the averages is to recenter each row of the fine camera resolution matrix.

Note that this change still does not produce a spectro-perfectionist resolution matrix. In other words, Rthis x model != observed.

julienguy commented 1 year ago

Code validated with the simulation of one delta emission line at 4000A : the resolution and output line rms are now identical. new-code

sbailey commented 1 year ago

It appears that specsim isn't running automated tests anymore (issue #129), and as a result it also appears that this PR broke the test_camera.py unit test.

It might be that the test itself needs to be updated. @p-slash please check the test and submit a new PR that either fixes the code or the test as needed:

cd (specsim_git_clone_directory)
pytest specsim/tests
...
====================================================== FAILURES =======================================================
__________________________________________________ test_downsampling __________________________________________________

    def test_downsampling():
        c = specsim.config.load_config('test')
        i = specsim.instrument.initialize(c)
        camera = i.cameras[0]

        # Use an intermediate dense matrix for downsampling.
        # This is the old implementation of get_output_resolution_matrix()
        # which uses too much memory.
        n = len(camera._output_wavelength)
        m = camera._downsampling
        i0 = camera.ccd_slice.start - camera.response_slice.start
        R1 = (camera._resolution_matrix[: n * m, i0 : i0 + n * m].toarray()
             .reshape(n, m, n, m).sum(axis=3).sum(axis=1) / float(m))

        # Use the new sparse implementation of get_output_resolution_matrix().
        R2 = camera.get_output_resolution_matrix()

>       assert np.allclose(R1, R2.toarray())
E       AssertionError: assert False
E        +  where False = <function allclose at 0x7f57ef063880>(array([[0.6386939 , 0.17828141, 0.00237037, ..., 0.        , 0.        ,\n        0.        ],\n       [0.17828141, 0.63...9 ,\n        0.17828141],\n       [0.        , 0.        , 0.        , ..., 0.0023727 , 0.17828141,\n        0.6386939 ]]), array([[7.21582136e-01, 1.38638623e-01, 5.69036350e-04, ...,\n        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n...\n       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n        5.71582129e-04, 1.38638623e-01, 7.21582136e-01]]))
E        +    where <function allclose at 0x7f57ef063880> = np.allclose
E        +    and   array([[7.21582136e-01, 1.38638623e-01, 5.69036350e-04, ...,\n        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n...\n       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n        5.71582129e-04, 1.38638623e-01, 7.21582136e-01]]) = <bound method spmatrix.toarray of <1500x1500 sparse matrix of type '<class 'numpy.float64'>'\n   with 7494 stored elements (5 diagonals) in DIAgonal format>>()
E        +      where <bound method spmatrix.toarray of <1500x1500 sparse matrix of type '<class 'numpy.float64'>'\n    with 7494 stored elements (5 diagonals) in DIAgonal format>> = <1500x1500 sparse matrix of type '<class 'numpy.float64'>'\n with 7494 stored elements (5 diagonals) in DIAgonal format>.toarray

specsim/tests/test_camera.py:40: AssertionError
p-slash commented 1 year ago

I've redesigned the test for the new resolution matrix in PR #132.