Bumped version.

962297f4 · Danilo Ferreira de Lima · 8f0ff399 · 962297f4 · 962297f4
Commit 962297f4 authored 1 year ago by Danilo Ferreira de Lima
--- a/notebook/Example offline analysis.ipynb
+++ b/notebook/Example offline analysis.ipynb
@@ -514,6 +514,15 @@
    "plt.ylim((None, 1.05))"
   ]
  },
+  {
+   "cell_type": "markdown",
+   "id": "d571f045-a877-42fd-9934-d9686bac4283",
+   "metadata": {},
+   "source": [
+    "### Resolution assessment per energy\n",
+    "\n"
+   ]
+  },
  {
   "cell_type": "markdown",
   "id": "ffd654f9-8175-45da-9ae1-09c184a7e520",

 %% Cell type:markdown id:59f50187-f73f-471b-b668-9126e6f48501 tags:

 # Learning high-resolution data from low-resolution

 This is an example notebook showing how to use the `pes_to_spec` infrastructure in this package.

 We start by importing some modules. The key module here is called `pes_to_spec`.

 %% Cell type:code id:d44af0b6-9c00-4e70-b49b-d74ed562e92f tags:

 ``` python
 import sys
 # add the pes_to_spec main directory
 # (change this depending on where you started the notebook if needed, or comment it out if you have done pip install in pes_to_spec)
 sys.path.append('..')

 # you meay need to do pip install matplotlib seaborn extra_data for this notebook, additionally
 # for this notebook the following packages are needed:
 # pip install "numpy>=1.21" "scipy>=1.6" "scikit-learn>=1.2.0" torch torchbnn  matplotlib seaborn extra_data
 ```

 %% Cell type:code id:da002d3e-c0da-419b-922b-0ab5c6deece8 tags:

 ``` python
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.gridspec import GridSpec
 import seaborn as sns

 import lmfit
 import scipy
 from extra_data import open_run, by_id
 from itertools import product
 from pes_to_spec.model import Model, matching_ids

 from typing import Any, Dict
 ```

 %% Cell type:markdown id:494a729c-dff4-4501-b828-fba2aaae5a23 tags:

 # Input data

 Read data from two runs. One shall be used for training the model. The second one is used for testing it.
 Note that the data in the training run must be large enough, compared to the number of model parameters.

 Only the SPEC, PES and XGM data is used for training, while only the PES and XGM data is needed for testing.
 However, more data is collected here to validate the results.

 %% Cell type:code id:4a301f2a-dedb-46e4-b096-fc9c6cf5b23a tags:

 ``` python
 run = open_run(proposal=900331, run=69)
 run_test = open_run(proposal=900331, run=70)

 # useful names to avoid repeating it all over the notebook, in case they ever change
 spec_name = "SA3_XTD10_SPECT/MDL/SPECTROMETER_SCS_NAVITAR:output"
 pes_name = "SA3_XTD10_PES/ADC/1:network"
 xgm_name = "SA3_XTD10_XGM/XGM/DOOCS:output"

 pres_name = "SA3_XTD10_PES/GAUGE/G30310F"
 volt_name = "SA3_XTD10_PES/MDL/DAQ_MPOD"

 # PES channels
 channels = [f"channel_{i}_{l}" for i, l in product([1, 3, 4], ["A", "B", "C", "D"])]

 def get_gas(run) -> str:
    """Get gas in chamber for logging."""
    gas_sources = [
                  "SA3_XTD10_PES/DCTRL/V30300S_NITROGEN",
                  "SA3_XTD10_PES/DCTRL/V30310S_NEON",
                  "SA3_XTD10_PES/DCTRL/V30320S_KRYPTON",
                  "SA3_XTD10_PES/DCTRL/V30330S_XENON",
              ]
    gas_active = list()
    for gas in gas_sources:
        # check if this gas source is interlocked
        if gas in run.all_sources and run[gas, "interlock.AActionState.value"].ndarray().sum() == 0:
            # it is not, so this gas was used
            gas_active += [gas.split("/")[-1].split("_")[-1]]
    gas = "_".join(gas_active)
    return gas

 def get_tids(run, need_spec:bool=True) -> np.ndarray:
    """Get which train IDs contain all necessary inputs for training."""
    spec_tid = run[spec_name, "data.trainId"].ndarray()
    pes_tid = run[pes_name, "digitizers.trainId"].ndarray()
    xgm_tid = run[xgm_name, "data.trainId"].ndarray()

    # match tids to be sure we have all inputs:
    tids = matching_ids(spec_tid, pes_tid, xgm_tid)
    return tids

 def get_data(run, tids) -> Dict[str, Any]:
    """Get all relevant data."""
    data = dict()
    data["int"] = run[xgm_name, "data.intensitySa3TD"].select_trains(by_id[tids]).ndarray()[:, 0][:, np.newaxis]
    data["pressure"] = run[pres_name, "value"].select_trains(by_id[tids]).ndarray()
    data["voltage"] = run[volt_name, "u212.value"].select_trains(by_id[tids]).ndarray()
    data["energy"] = run[spec_name, "data.photonEnergy"].select_trains(by_id[tids]).ndarray()
    data["spec"] = run[spec_name, "data.intensityDistribution"].select_trains(by_id[tids]).ndarray()
    data["pes"] = {ch: run[pes_name,
                           f"digitizers.{ch}.raw.samples"].select_trains(by_id[tids]).ndarray()
                    for ch in channels}
    data["gas"] = get_gas(run)
    return data
 ```

 %% Cell type:code id:210c0550-1abb-43a0-99a5-7c35d2766be0 tags:

 ``` python

 # get the matched train IDs
 tids = get_tids(run)

 # we don't need the spec for testing in reality,
 # but it is nice to plot it in the test run too,
 # to check that this works during validation
 test_tids = get_tids(run_test, need_spec=True)

 # get the data
 data = get_data(run, tids)
 data_test = get_data(run_test, test_tids)
 ```

 %% Cell type:markdown id:017865a1-057f-48c7-8bef-a6e40490de2c tags:

 Now the `data` and `data_test` dictionaries contain the necessary information about the training and test runs.
 The code above also selected only entries with train IDs on which at least SPEC, PES and XGM were present.

 Note that for training, it is assumed that only one pulse is present. For testing there is no such requirement.

 First output some general information about the conditions of the measurement device.

 %% Cell type:code id:956105a6-d37e-453c-bfeb-2b1c876ee3f2 tags:

 ``` python
 print(f"Gas in training: {data['gas']}")
 print(f"Gas in testing: {data_test['gas']}")
 ```

 %% Output

    Gas in training: NEON
    Gas in testing: NEON

 %% Cell type:code id:4654f205-edc6-45f7-97bd-0d088c38edb0 tags:

 ``` python
 print(f"Voltage in training: {np.mean(data['voltage']):0.2f} +/- {np.std(data['voltage']):0.2f}")
 print(f"Voltage in testing: {np.mean(data_test['voltage']):0.2f} +/- {np.std(data_test['voltage']):0.2f}")
 ```

 %% Output

    Voltage in training: -116.00 +/- 0.01
    Voltage in testing: -116.00 +/- 0.01

 %% Cell type:code id:fa662544-3caa-4404-bb61-fa41add82642 tags:

 ``` python
 print(f"Pressure in training: {np.mean(data['pressure']):0.2e} +/- {np.std(data['pressure']):0.2e}")
 print(f"Pressure in testing: {np.mean(data_test['pressure']):0.2e} +/- {np.std(data_test['pressure']):0.2e}")
 ```

 %% Output

    Pressure in training: 1.29e-06 +/- 3.95e-08
    Pressure in testing: 1.29e-06 +/- 3.91e-08

 %% Cell type:markdown id:d9b62e3a-aff9-4436-905d-ea02c11aecf6 tags:

 ## Using the virtual spectrometer

 %% Cell type:markdown id:5962e483-60da-4c70-bb09-dce5fc9745e0 tags:

 Now we will actually train the model. We do that by creating a `Model` object (from `pes_to_spec`) and calling the `fit` function.
 The `fit` function requires the PES intensity, the SPEC intensity, the energy axis from SPEC (stored as a reference only), as well as the energy measured in the XGM (which has better resolution than the integral of the PES).

 %% Cell type:code id:a0adb57b-7496-4781-9511-ac2a8d05658d tags:

 ``` python
 # this is the main object holding all
 # information needed for training and prediction
 # the default parameters should be sufficient in most times
 model = Model(channels=channels,
              high_res_sigma=0.0,
             )

 # this trains the model
 # the first parameter is expected to be a dictionary with the channel name as a key
 model.fit(data['pes'],
          data['spec'],
          data['energy'],
          pulse_energy=data['int'])

 # save it for later usage:
 model.save("model.joblib")

 # load a model (you can start from here if working on an existing model)
 model = Model.load("model.joblib")

 # and use it to map a low-resolution spectrum to a high-resolution one
 # as before, the low_resolution_raw_data refers to a dictionary mapping the channel name
 # in the format "channel_[1-4]_[A-D]" to the 2D numpy array with shape (number_of_train_IDs, features)
 # all names and shapes must match the format in training, except for the number_of_train_IDs, which may vary
 pred = model.predict(data['pes'], pulse_energy=data['int'])
 ```

 %% Output

    Checking data quality in high-resolution data.
    Selected 7165 of 7165 samples
    Fitting PCA on low-resolution data.
    Using 1000 comp. for PES PCA (asked for 1000, out of 7201, in 7165 samples).
    Fitting PCA on high-resolution data.
    Fitting outlier detection
    Fitting model.
    Calculate PCA unc. on high-resolution data.
    Calculate transfer function
    Resolution: 0.9825132876894322
    Calculate PCA on channel_1_A
    Calculate PCA on channel_1_B
    Calculate PCA on channel_1_C
    Calculate PCA on channel_1_D
    Calculate PCA on channel_3_A
    Calculate PCA on channel_3_B
    Calculate PCA on channel_3_C
    Calculate PCA on channel_3_D
    Calculate PCA on channel_4_A
    Calculate PCA on channel_4_B
    Calculate PCA on channel_4_C
    Calculate PCA on channel_4_D
    End of fit.

 %% Cell type:markdown id:e0286ae3-1a59-468f-ae40-c3ed94b7b301 tags:

 Now we can try it in the test dataset:

 %% Cell type:markdown id:ffc06362-3479-4cb9-b102-b438a83d2950 tags:

 We can predict it in the training data itself, but this is a bit biased, since we used the same information to fit the model.

 %% Cell type:code id:917156f3-9476-48e0-9121-5f75f185045f tags:

 ``` python
 pred = model.predict(data_test['pes'], pulse_energy=data_test['int'])

 # add the references in this array in the same array format, so we can plot them later
 pred["energy"] = model.get_energy_values()

 pred['spec'] = data_test['spec'][:, np.newaxis, :]
 ```

 %% Cell type:markdown id:77866435-1cb6-40ac-a9a5-8f2ae37017b7 tags:

 Let's try to predict in the independent run in the test dataset. The performance of the model varies a lot if the beam intensity is very different from the training one. To ensure we take a train ID to visualize that is relatively high intensity, we sort the train IDs by XGM intensity and then choose the highest intensity one.
 One could try other train IDs.

 For train IDs with close to zero beam intensity, there is a relatively larger error, since the training data did not contain any of those samples and the signal-to-noise ratio is relatively high.

 %% Cell type:code id:ed62606a-4ea7-4e0a-8b61-73e682cacf04 tags:

 ``` python
 # choose train ID of the test dataset by XGM intensity
 test_intensity = np.argsort(data_test['int'][:,0])
 example_tid = test_intensity[-1]
 ```

 %% Cell type:markdown id:f931e9e0-84a7-4e4f-bfad-588bfe77267c tags:

 Now we can actually plot it.

 %% Cell type:code id:fd42984c-554c-4c69-bf8a-119eeb0cca62 tags:

 ``` python
 def plot(data):
    """Plot prediction and expectation."""
    fig = plt.figure(figsize=(12, 8))
    gs = GridSpec(1, 1)
    ax = fig.add_subplot(gs[0, 0])
    ax.plot(data["energy"], data["spec"], c='b', lw=3, label="High-res. measurement")
    ax.plot(data["energy"], data["expected"], c='r', ls='--', lw=3, label="High-res. prediction")
    ax.fill_between(data["energy"], data["expected"] - data["total_unc"], data["expected"] + data["total_unc"], facecolor='gold', alpha=0.5, label="68% unc.")
    ax.legend(frameon=False, borderaxespad=0, loc='upper left')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    Y = np.amax(data["spec"])
    ax.set(
            xlabel="Photon energy [eV]",
            ylabel="Intensity",
            ylim=(0, 1.3*Y))
    plt.show()
 ```

 %% Cell type:code id:bbbf77b5-f914-4b47-8ab6-fd3a89d0f983 tags:

 ``` python
 # select the correct train ID for the data to plot
 # except for the energy axis, which is always the same
 plot({k: v[example_tid, 0, :] if k != "energy" else v
      for k, v in pred.items()
      if k in ["expected", "total_unc", "spec", "energy"]})
 ```

 %% Output



 %% Cell type:markdown id:6a9c996e-1462-45a3-a9af-7ccbcc37ad69 tags:

 ## Further studies and validation

 The next items show methods to estimate the resolution and improve the results.

 They are not done by default, as they make further assumptions, which cannot always be applied. The code for such analyses is shown here, for reference, but one should be careful about how they are used.

 %% Cell type:markdown id:eca3a06a-d613-4206-b128-6d81031de1d1 tags:

 ### Resolution assessment using the autocorrelation

 We establish the resolution of the virtual spectrometer using the autocorrelation function, which estimates which level of detail can be observed in the test dataset.

 The autocorrelation function cannot assess which effect are physically relevant and which are simply noise. Therefore this method can only provide a rough estimate of the resolution. It is not expected to be very precise, but it can be used for a quick assessment.

 %% Cell type:code id:1491550c-6940-425e-a557-f2cd381d287b tags:

 ``` python
 def fwhm(x: np.ndarray, y: np.ndarray) -> float:
    """Return the full width at half maximum of x."""
    # half maximum
    half_max = np.amax(y)*0.5
    # signum(y - half_max) is zero before and after the half maximum,
    # and it is 1 in the range above the half maximum
    # The difference will be +/- 1 only at the transitions
    d = np.diff(np.sign(y - half_max))
    left_idx = np.where(d > 0)[0][0]
    right_idx = np.where(d < 0)[-1][-1]
    return x[right_idx] - x[left_idx]

 def autocorrelation(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Given the energy axis in x and the intensity in y, calculate the auto-correlation function."""
    mean_y = np.mean(y, keepdims=True, axis=0)
    e = x - np.mean(x)
    Rxx = np.mean(np.fft.fftshift(np.fft.ifft(np.absolute(np.fft.fft(y - mean_y))**2), axes=(-1,)), axis=(0,1))
    Rxx /= np.amax(Rxx)
    return Rxx
 ```

 %% Cell type:code id:45fc52a4-5716-42b9-adbd-a307d16c0c34 tags:

 ``` python
 fig = plt.figure(figsize=(10, 8))
 R = dict()
 res = dict()
 for instr, title in {"expected": "Virtual spectrometer",
                     "spec":"Grating spectometer",
                     #"pes": "PES",
                    }.items():
    e = pred["energy"] - np.mean(pred["energy"])
    R[instr] = autocorrelation(pred["energy"], pred[instr])
    res[instr] = fwhm(e, R[instr])
    plt.plot(e, R[instr], lw=2, label=f"{title} (FWHM = {res[instr]:.2f} eV)")

 plt.legend(frameon=False)
 plt.xlabel("Energy [eV]")
 plt.ylabel("Autocorrelation")
 plt.xlim((-3, 3))
 plt.ylim((None, 1.05))
 ```

 %% Output

    (-0.17591992514901664, 1.05)



+%% Cell type:markdown id:d571f045-a877-42fd-9934-d9686bac4283 tags:
+
+### Resolution assessment per energy
+
+
 %% Cell type:markdown id:ffd654f9-8175-45da-9ae1-09c184a7e520 tags:

 ### Resolution assessment using deconvolution

 Here we attempt to establish the resolution of the virtual spectrometer using a deconvolution-based method. The idea here is that the virtual spectrometer can be seen as a *linear* device that somehow *worsens* the resolution of the grating spectrometer. Within the context of linear systems theory any such device can be modelled mathematically as a block that applies a convolution between a function $g$ and the grating spectrometer data.

 That is, if the grating spectrometer data is $y$ and the virtual spectrometer result is $\hat{y}$, then we assume that there is a function $g$ such that:

 $\hat{y} = y \ast g + \epsilon$,

 where $\epsilon$ is zero-mean Gaussian noise.

 Under such an approach, one can calculate the function $g$ exactly, by performing a deconvolution between $\hat{y}$ and $y$.

 %% Cell type:code id:7ed071e5-4f60-4195-830a-73ab8e5c2577 tags:

 ``` python
 def deconv(y: np.ndarray, yhat: np.ndarray) -> np.ndarray:
    """Given the grating spectrometer data and the virtual spectrometer data,
    calculate the deconvolution between them.
    """
    # subtract the mean spectra to remove the FEL bandwidth
    yhat_s = yhat - np.mean(yhat, keepdims=True, axis=(0, 1))
    y_s = y  - np.mean(y, keepdims=True, axis=(0, 1))
    # Fourier transforms
    Yhat = np.fft.fft(yhat_s)
    Y = np.fft.fft(y_s)
    # spectral power of the assumed "true" signal (the grating spectrometer data)
    Syy = np.mean(np.absolute(Y)**2, axis=(0, 1))
    Syh = np.mean(Y*np.conj(Yhat), axis=(0, 1))
    # approximate transfer function as the ratio of power spectrum densities
    H = Syh/Syy
    return np.fft.fftshift(np.fft.ifft(H))
 ```

 %% Cell type:code id:a1c5137f-fe6b-4930-aff5-90026ad5f3c3 tags:

 ``` python
 # centered energy axis
 e = pred["energy"] - np.mean(pred["energy"])
 # impulse response
 g = deconv(pred["spec"], pred["expected"])
 ```

 %% Cell type:code id:75d3ddc3-a4b6-4869-bc1c-f08320089845 tags:

 ``` python
 def fit_gaussian(x: np.ndarray, y: np.ndarray):
    """Fit Gaussian."""
    def gaussian(x, amp, cen, wid):
        return amp * np.exp(-0.5 * (x-cen)**2 / (wid**2))
    gmodel = lmfit.Model(gaussian)
    result = gmodel.fit(y, x=x, cen=0.0, amp=1.0, wid=1.0)
    return result.best_values["wid"]*2.355, result
 ```

 %% Cell type:code id:f9bbd13c-d972-4af1-a6c1-0fe12509317a tags:

 ``` python
 width, result = fit_gaussian(e, np.absolute(g))
 ```

 %% Cell type:code id:26641ed6-47cd-418d-ab3c-e63c83962387 tags:

 ``` python
 plt.figure(figsize=(10, 8))
 plt.plot(e, np.absolute(g), lw=2, label="Virtual spectrometer")
 plt.plot(e, result.best_fit, lw=2, label=f"Gaussian fit, FWHM = {width:.2f} eV")
 plt.xlim(-3, 3)
 plt.xlabel("Energy [eV]")
 plt.ylabel("Impulse response [a.u.]")
 plt.legend(frameon=False)
 ```

 %% Output

    <matplotlib.legend.Legend at 0x2b5378c4e700>



 %% Cell type:markdown id:00a1bdb9-b52a-4f8b-8c03-30407f486595 tags:

 Note that this response function does *not* tell us the resolution of the virtual spectrometer. It tells us how we can smear the grating spectrometer data to transform that data into the virtual spectrometer. That is, this is how much worse we do with the virtual spectrometer, relative to the grating spectrometer.

 As a result, if we approximate the response functions with Gaussians and assume that the previous autocorrelation function gives us an estimate of the grating spectrometer resolution, we can guess the total resolution as:

 $\sigma_{total} = \sqrt{\sigma_{grating}^2 + \sigma_{VIRT}^2}$

 The same relation is applies for the FWHM. This relation assumes independence between the two systems and assumes we can approximate the response functions as Gaussians.

 %% Cell type:code id:2f466226-fbad-4f64-86d3-9b5d0685678b tags:

 ``` python
 total_resolution = np.sqrt(width**2 + res["spec"]**2)
 ```

 %% Cell type:code id:051d329e-3527-4ce1-abfe-6da441ce2e4b tags:

 ``` python
 total_resolution
 ```

 %% Output

    1.0442815003593549

 %% Cell type:markdown id:2c85f6fd-c2a5-419f-ae61-8fe77289cf7d tags:

 The previously obtained resolution of the virtual spectrometer using the autocorrelation method was:

 %% Cell type:code id:e8d3b2f8-09ba-4e9f-9b54-b48867dadc89 tags:

 ``` python
 res["expected"]
 ```

 %% Output

    0.9825132876894713

 %% Cell type:markdown id:efac989d-34e1-4378-abe1-b215bf58c05c tags:

 Notice, however, that the response function is not Gaussian and therefore, one could use the full function. to actually simulate the virtual spectrometer.

 Furthermore, this ignores the uncertainty effect, which could be seen as an extra noise level added on top of the virtual spectrometer.

 %% Cell type:markdown id:dd1f6722-472f-488d-a0eb-e362732fd364 tags:

 ### Validation: compare grating spectrometer and simulated virtual spectrometer

 To check that the resolution estimate is correct, we take an example grating spectrometer pulse and smear it by the impulse response function $g$ above. If it is correct, we should get a similar result as the virtual spectrometer itself.

 %% Cell type:code id:36e25952-b295-43d1-bffc-d62dc1cd0a95 tags:

 ``` python
 # smearing
 sigma = width/2.355 # convert FWHM to sigma
 # create a Gaussian with the requested width
 g_simple = np.exp(-0.5 * (pred["energy"] - np.mean(pred["energy"]))**2/(sigma**2))
 g_simple  /= np.sum(g_simple)
 # smear the grating spectrometer data
 y_simul = scipy.signal.fftconvolve(pred["spec"], g_simple*np.ones_like(pred["spec"]), mode="same", axes=-1)
 ```

 %% Cell type:code id:e25f2828-fb32-4747-b7ea-24752c390625 tags:

 ``` python
 fig = plt.figure(figsize=(10, 8))
 plt.plot(pred["energy"], y_simul[example_tid,0], c='b', lw=3, label="Smeared grating spec.")
 plt.plot(pred["energy"], pred["expected"][example_tid,0], c='r', ls='--', lw=3, label="Virtual spectrometer")
 plt.fill_between(pred["energy"],
                 pred["expected"][example_tid, 0] - pred["total_unc"][example_tid,0],
                 pred["expected"][example_tid,0] + pred["total_unc"][example_tid,0],
                 facecolor='gold', alpha=0.5, label="68% unc.")
 plt.xlabel("Energy [eV]")
 plt.ylabel("Intensity [a.u.]")
 plt.legend(frameon=False)
 ```

 %% Output

    <matplotlib.legend.Legend at 0x2b53789a2ca0>



 %% Cell type:markdown id:db3ddff6-dad5-40f5-af6c-680ca2657a24 tags:

 ## Improve the resolution further: Wiener deconvolution

 If we know the impulse response of the virtual spectrometer, we can undo that effect. This assumes however, that the resolution function is very accurate. This may not be true, as approximations are made previously (such as assuming the same resolution for all energies and linearity).

 Given the limitation created by the uncertainty, this is often not very reliable.

 %% Cell type:code id:cb27ad10-5cb0-4f77-a6d6-8a1a41b6853f tags:

 ``` python
 dec = model.deconvolve(pred["expected"])
 ```

 %% Cell type:code id:a52d01fb-878f-4f3b-940f-c47013df24ca tags:

 ``` python
 fig = plt.figure(figsize=(10, 8))
 plt.plot(pred["energy"], pred["spec"][example_tid,0], c='b', lw=3, label="Grating spec.")
 plt.plot(pred["energy"], dec[example_tid,0], c='r', ls='--', lw=3, label="Virtual spectrometer (deconvolved)")
 plt.fill_between(pred["energy"],
                 dec[example_tid, 0] - pred["total_unc"][example_tid,0],
                 dec[example_tid,0] + pred["total_unc"][example_tid,0],
                 facecolor='gold', alpha=0.5, label="68% unc.")
 plt.xlabel("Energy [eV]")
 plt.ylabel("Intensity [a.u.]")
 plt.legend(frameon=False)
 ```

 %% Output

    <matplotlib.legend.Legend at 0x2b56166ab340>



 %% Cell type:code id:76eaa0ed-d4c5-426d-9b10-c58b14ffb059 tags:

 ``` python
 ```

--- a/pes_to_spec/__init__.py
+++ b/pes_to_spec/__init__.py
@@ -2,4 +2,4 @@
 Estimate high-resolution photon spectrometer data from low-resolution non-invasive measurements.
 """

-VERSION = "0.3.2"
+VERSION = "0.3.3"