simulators with alternate wavelength grids

[sbailey]

Test driven development: This PR provides a test that demonstrates a bug in calculating the pixel_edges vs. sim_edges for the wavelength binning. My initial attempt at refactoring int vs. round vs. integer division just swapped what succeeded vs. failed, so I'm submitting this report with a test to get more eyes on it.

The test currently fails, but presumably we will find a fix and add that before merging this.

============================================================================ FAILURES ============================================================================
______________________________________________________________________ test_alt_wavelengths ______________________________________________________________________

    def test_alt_wavelengths():
        import yaml
        import pkg_resources
        configfile = pkg_resources.resource_filename(
            'specsim', 'data/config/test.yaml')
        configdata = yaml.load(open(configfile))

        # Test input 0.1 -> output 1.2 Angstrom wavelength grid
        configdata['wavelength_grid']['step'] = 0.1
        configdata['instrument']['cameras']['r']['constants']['output_pixel_size'] = '1.2 Angstrom'
        config = specsim.config.Configuration(configdata)
        sim = Simulator(config)

        # Test input 0.1 -> output 0.6 Angstrom wavelength grid
        configdata['wavelength_grid']['step'] = 0.1
        configdata['instrument']['cameras']['r']['constants']['output_pixel_size'] = '0.6 Angstrom'
        config = specsim.config.Configuration(configdata)
        sim = Simulator(config)

        # Test input 0.1 -> output 0.4 Angstrom wavelength grid
        configdata['wavelength_grid']['step'] = 0.1
        configdata['instrument']['cameras']['r']['constants']['output_pixel_size'] = '0.4 Angstrom'
        config = specsim.config.Configuration(configdata)
>       sim = Simulator(config)

specsim/tests/test_simulator.py:46: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
specsim/simulator.py:66: in __init__
    self.instrument = specsim.instrument.initialize(config)
specsim/instrument.py:606: in initialize
    constants['num_sigmas_clip'], constants['output_pixel_size']))
specsim/camera.py:236: in __init__
    pixel_edges, sim_edges, rtol=0., atol=1e-6 * wavelength_step):
/Users/sbailey/anaconda/envs/desi/lib/python3.5/site-packages/numpy/core/numeric.py:2372: in allclose
    res = all(isclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan))
/Users/sbailey/anaconda/envs/desi/lib/python3.5/site-packages/numpy/core/numeric.py:2453: in isclose
    return within_tol(x, y, atol, rtol)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

x = array([ 5599.95,  5600.35,  5600.75, ...,  7398.75,  7399.15,  7399.55]), y = array([ 5599.95,  5600.35,  5600.75, ...,  7399.15,  7399.55,  7399.95])
atol = 9.9999999999909045e-08, rtol = 0.0

    def within_tol(x, y, atol, rtol):
        with errstate(invalid='ignore'):
>           result = less_equal(abs(x-y), atol + rtol * abs(y))
E           ValueError: operands could not be broadcast together with shapes (4500,) (4501,)

/Users/sbailey/anaconda/envs/desi/lib/python3.5/site-packages/numpy/core/numeric.py:2436: ValueError
======================================================== 1 failed, 71 passed, 9 skipped in 22.91 seconds =========================================================

[sbailey]

Tests pass on my laptop, though still waiting for Travis. Assuming they pass on Travis too, this is ready for review.

[dkirkby]

There are two separate issues here:

• how should users change the default wavelength grid?
• reasonable user-defined grids do not work as expected.

It looks like your changes have fixed the second issue, but I would still like to understand the first issue, and avoid hot wiring the YAML in memory if possible ;-)

What is the use case for changing the default grid? The grid needs to be specified before the various objects that make up a simulator are initialized, since they all interpolate to that grid. But I'm guessing you don't need to make changes after the initialization phase?

[sbailey]

I'd be happy to have a different way to update the default wavelength grid, but I'd still like to get the wavelength grid fix in soon since otherwise I'm stuck working off my branch in the meantime for mocks -> spectra development.

Use case for changing default grid: select_mock_targets writes out truth spectra for all targets on some grid. I need to read that in and then create a simulator that matches that grid. I don't need to change the grid once the simulator is created so changing the default grid at the level of the Configuration object before making a Simulator is fine.

Other notes:

• I don't want to create a custom desi.yaml file that exactly matches the mock grid because that tightly couples the mock generation to the simulations of those mocks. Introspecting the mocks to discover the grid and make a matching simulator is more flexible.
• In general we should be able to update any config parameter in memory without having to go via the file. i.e. let's not treat updating the wavelength grid as a special case. I originally tried manipulating the Configuration object itself after creating it from the yaml file but couldn't work out the recipe, then I resorted to editing the yaml structure in memory ahead of time like the tests do, and then I found that my requisite wavelength grid (0.2 A input and output) didn't work anyway, thus this PR.

[dkirkby]

I suggest the following steps to change the wavelength grid:

desi_config = specsim.config.load_config('desi')
desi_config.wavelength_grid.min = 4000.
desi_config.wavelength_grid.max = 5000.
desi_config.wavelength_grid.step = 100.
desi_config.update()

You can check that this does the right thing using:

print desi_config.wavelength
[ 4000.  4100.  4200.  4300.  4400.  4500.  4600.  4700.  4800.  4900.   5000.] Angstrom

Now build a simulator with the updated config:

desi = specsim.simulator.Simulator(desi_config)

[dkirkby]

The general philosophy of specsim is that the simulation grid should be selected to balance accuracy and speed in the simulation and all simulation inputs should be specified on their own natural grid. The config initialization mechanism internally resamples all inputs to the simulation grid, but I didn't implement that yet for the spectra (or other functions of wavelength) passed directly to the updated simulate() method (which bypasses the config mechanism).

[dkirkby]

Can you update the unit tests to use the method above for changing the grid step size and then lets merge this now and revisit resampling of simulate() inputs later.

[sbailey]

The unit tests include cases of also changing the output grid size so I can't update them to this method without knowing the proper recipe for that too. That's also what I need for my simulation use case.

Here are some examples of getting tripped up trying to change the output grid size:

desi_config = specsim.config.load_config('desi')
desi_config.wavelength_grid.min = 3500.
desi_config.wavelength_grid.max = 7000.
desi_config.wavelength_grid.step = 100.

desi_config.instrument.cameras.b.constants.output_pixel_size = 100
desi_config.instrument.cameras.r.constants.output_pixel_size = 100
desi_config.instrument.cameras.z.constants.output_pixel_size = 100

desi_config.update()
specsim.simulator.Simulator(desi_config)

Causes "UnitConversionError: '' (dimensionless) and 'Angstrom' (length) are not convertible" at the time of creating the Simulator.

Trying with units:

desi_config = specsim.config.load_config('desi')
desi_config.wavelength_grid.min = 3500.
desi_config.wavelength_grid.max = 7000.
desi_config.wavelength_grid.step = 100.

desi_config.instrument.cameras.b.constants.output_pixel_size = 100 * u.Angstrom
desi_config.instrument.cameras.r.constants.output_pixel_size = 100 * u.Angstrom
desi_config.instrument.cameras.z.constants.output_pixel_size = 100 * u.Angstrom

desi_config.update()
specsim.simulator.Simulator(desi_config)

Causes "TypeError: Only dimensionless scalar quantities can be converted to Python scalars"

And finally I noticed that before I try to alter it desi_config.instrument.cameras.b.constants.output_pixel_size is still a string, so let's try that:

desi_config = specsim.config.load_config('desi')
desi_config.wavelength_grid.min = 3500.
desi_config.wavelength_grid.max = 7000.
desi_config.wavelength_grid.step = 100.

desi_config.instrument.cameras.b.constants.output_pixel_size = "100 Angstrom"
desi_config.instrument.cameras.r.constants.output_pixel_size = "100 Angstrom"
desi_config.instrument.cameras.z.constants.output_pixel_size = "100 Angstrom"

desi_config.update()
specsim.simulator.Simulator(desi_config)

Still doesn't work: "IndexError: index 0 is out of bounds for axis 0 with size 0"

That's when I resorted to hotwiring the input yaml in memory... :)

[dkirkby]

Your last attempt looks fine to me (since config constants are always strings, as you noticed), but fails (with a helpful exception message) because your min/max do not cover the range of non-zero camera throughputs in $DESIMODEL. However, it works for me with the change:

desi_config.wavelength_grid.max = 9900.

[sbailey]

I was closer than I realized. I updated the tests to use that method for altering input and output wavelength grids. I think this addresses both the initial fix for the wavelength binning and tests it via the recommended method of altering the wavelength grids of the configuration without hotwiring yaml structures on the fly.

[dkirkby]

Merging now, thanks!