Add SpectrumBeam as derived class from Beam object

[phyy-nx]

SpectrumBeam has these methods:
set_spectrum(flex.double energies, flex.double weights)
get_spectrum_energies()
get_spectrum_weights()
get_weighted_wavelength()

Two flex arrays, energies (eV) and weights (dimensionless), represent the spectrum as a scatter plot. This follows the conventions set up in nexusformat/definitions#717. The new function get_weighted_wavelength uses a simple weighted mean to determine a wavelength.

It is ok if the SpectrumBeam has a wavelength != get_weighted_wavelength. This represents a calibrated wavelength determined from processing the spectrum.

is_similar_to tests get_weighted_wavelength between the two beams.

Includes tests.

Marking as a draft until #151 is merged, then will rebase on master and mark as ready for review. This is ready for review now though.

[phyy-nx]

I'm not quite sure what to do with is_similar_to when a Beam is compared to a SpectrumBeam. As noted, when two SpectrumBeams are compared to eachother, get_weighted_wavelength is used as a test. That's fine. However, here is the behavior now between Beam and SpectrumBeam (paraphrased):

a = Beam()
b = SpectrumBeam()
a.is_similar_to(b) # works because of polymorphim
True
b.is_similar_to(a) # uses dynamic cast
Error (a has no spectra)

Certainly it should be commutative, but should it be true, false, or throw an exception? I lean towards false.

[graeme-winter]

I'm not quite sure what to do with is_similar_to when a Beam is compared to a SpectrumBeam. As noted, when two SpectrumBeams are compared to eachother, get_weighted_wavelength is used as a test. That's fine. However, here is the behavior now between Beam and SpectrumBeam (paraphrased):

a = Beam()
b = SpectrumBeam()
a.is_similar_to(b) # works because of polymorphim
True
b.is_similar_to(a) # uses dynamic cast
Error (a has no spectra)

Certainly it should be commutative, but should it be true, false, or throw an exception? I lean towards false.

I would say they are not similar so agree a != b and b != a 🙂

[dagewa]

The way we went for the crystal recalculated cell was to only perform that part of the is_similar_to check if both models have the extended attribute. In that case we ended up useing boost::optional rather than inheritance, but perhaps the same logic applies here?

[phyy-nx]

Ok, with derived classes you have to use type_info to get comparators between derived class to work properly, which is outside of my scope for the moment. For now, I've made it so both variants of is_similar_to return true if the base class's is_similar_to returns true. So at least is it commutative.

Also checked in serialization fix. Good for more review. Thanks.

[dermen]

Should this functionality be merged with flex_Beam ?

from dxtbx_model_ext import flex_Beam
from dxtbx.model.beam import BeamFactory

beams = flex_Beam()
for wave,flux in [(1,1e12),(1.1,2e12),(1,1.5e12)]:
    beam = BeamFactory.simple(wave)
    beam.set_flux(flux)
    beams.append(beam)

print( [b.get_wavelength() for b in beams])
# [1.0, 1.1, 1.0]
print( [b.get_flux() for b in beams])
# [1000000000000.0, 2000000000000.0, 1500000000000.0]

[dagewa]

Making this a derived class results in a lot of boilerplate, where you have to duplicate the various constructors etc. of Beam. If Beam changes in the future the changes to constructors would have to be written in two places, or more if in future there are additional derived classes written.

There are alternatives. The simplest would be to put the new members and methods into Beam directly, leaving the data unset in most situations so it remains cheap. Another would be to use composition, whereby SpectrumBeam takes an already-constructed Beam, avoiding the need to rewrite all the constructors for this.

I don't have a strong opinion on the right approach, but I wanted to open this up for discussion.

[dagewa]

Perhaps put a comment here to note that this is the eV per Å conversion factor

[dagewa]

It would be better to append to the existing output for Beam so that if it changes, the output for SpectrumBeam will be updated automatically

[phyy-nx]

@dermen I wasn't aware of flex_Beam but I can see how it's useful. In this case I want to connect a spectrum with a specific experiment and serialize that through the experiment list. I don't think I could add a flex_Beam as the beam object of the Experiment without modifying flex_Beam to be polymorphic with Beam.

@dagewa we talked about your comment in group meeting today and there was general agreement. Polymorphism is rad, but lots of boilerplate does make maintenance harder. I'll look to see if the boilerplate can be minimized, but if not I think I'll rewrite this as members of Beam instead of using derived classes.

[dermen]

Having been working with these for several days now Im firmly opposed to the serialization of the spectra data (as we discussed), but I think all agree. Also I want to mention something so its in mind when developing the new code - the energies array should not change for each SpectrumBeam, just the weights. This is more to do with the storage of the raw spectra data however. So get_spectrum_energies() should always point to the same array. This is in line with everything we know about recording spectra at facilities, where the spectrometer doesnt actually change its number of pixels (which corresponds to the number of points along an energy axis). The thing that can and will change is the conversion of that energy axis in pixels to some units of keV, which is represented by a set of a few coefficients (at least 2 for a linear conversion). This is also good because it halves the raw data. Another way to think about it

# assume some arbitrary spectrum weights (doesnt matter for this)
num_energy = 10
weights = np.random.random(num_energy)

# assume the spectrum energies are
energies = np.linspace(9000, 9100, num_energy)
print(energies)

# array([9000.        , 9011.11111111, 9022.22222222, 9033.33333333,
       9044.44444444, 9055.55555556, 9066.66666667, 9077.77777778,
       9088.88888889, 9100.        ])

# You never really need to store the energies array, you just need two coefficients

energies2 = (9100 - 9000) / 9 * np.arange(num_energy) + 9000
assert np.allclose(energies2, energies)
#so you can just store
coeffs = (9100 - 9000) / 9,  9000

So I think best design would be if SpectrumBeam has members "weights" , "coefficients" and "num_energies" and then the method get_spectrum_energies applies the mapping to return the energies.

[graeme-winter]

You mention in #158 that you want to merge that one and this one, but this one is draft and @dermen raises some interesting & valid points.

I would comment here my brief experience of looking at X-FEL spectrometer data suggests that they are not always binned the same way, as they are themselves models of the pixel counts on the spectrometer - though I think serialising as a rectangular grid would be good.

I was also wondering on what scale the values are? e.g. should the sum be 1 or some arbitrary / unknown amount?

[phyy-nx]

Thanks, several questions here. First, this pull request will remain as a draft as I like the implementation in #158 better. Once #158 is merged, I'll close this one.

Regarding @dermen's comments. First, I agree with @graeme-winter that the binning of the xfel energy axis is not always a simple linear computation, so coefficients wouldn't work, much as I like the idea. Regarding storing a 2D spectra, the raw un-flattened spectrometer image, that's not really part of all this. For now, I think we are only interested in provided a 1D spectra.

(I'll put serialization/file sizes discussion in #158).

[graeme-winter]

It was two dimensional spectra in the sense of a one dimensional spectrum for each image I was referring to. That’s two dimensional data 🙂

[phyy-nx]

Oh gotcha. Sorry. For the scale, I think I'd prefer just using whatever the beamline recorded? Any normalization could be done by downstream applications.

[phyy-nx]

Superceded by #158