Allow frequency scaling (#5)

* Add initial frequency scaling

* Update tests for new parameter

* Refactor tests

- Distinguish between 'opts' (used as model inputs) and
 'pars' (which are used for scaling/fitting)
- Use parametrize rather than explicit loops for run_opts
- Use parametrize for materials rather than putting them in opts
- Make test options a list rather than a generator (opts must be
  iterated over multiple times when using parametrize)
- Simplify logic for removing None values from opts

* Add test for scaling parameters

* Add depreaction warnings

* Use kwargs rather than list of pars in horace_disp

* Update changelog

* Test positional and keyword argument combinations

* Multiple atol by scaling

CHANGELOG.rst

@@ -1,6 +1,16 @@
 `Unreleased <https://github.com/pace-neutrons/euphonic_sqw_models/compare/v0.4.0...HEAD>`_
 ----------
 
+- Improvements:
+
+  - There is a new ``frequency_scale=1.0`` argument to ``horace_disp`` which
+    allows the output frequencies to be scaled
+
+- Breaking changes:
+
+  - The ``pars=[]`` argument to ``horace_disp`` has been changed to
+    ``intensity_scale=1.0``
+
 `v0.4.0 <https://github.com/pace-neutrons/euphonic_sqw_models/compare/v0.3.0...0.4.0>`_
 ------
 

euphonic_sqw_models/euphonic_wrapper.py

@@ -16,10 +16,9 @@ class CoherentCrystal(object):
 
     Methods
     -------
-    w, sf = horace_disp(qh, qk, ql, pars=[])
+    w, sf = horace_disp(qh, qk, ql, intensity_scale=1.0, frequency_scale=1.0)
         Calculates the phonon dispersion surfaces `w` and structure factor `sf` at the specified
-        (qh, qk, ql) points. The optional parameter `pars=[scale_factor]` is just a single scale
-        factor to multiply the calculated structure factor by.
+        (qh, qk, ql) points, with optional intensity and frequency scaling factors
 
     Attributes
     ----------
@@ -110,17 +109,19 @@ def _calculate_sf(self, qpts):
             sf = np.hstack((sf, neg_sf))
         return w, sf
 
-    def horace_disp(self, qh, qk, ql, pars=[], *args, **kwargs):
+    def horace_disp(self, qh, qk, ql, intensity_scale=1.0, frequency_scale=1.0, *args, **kwargs):
         """
         Calculates the phonon dispersion surface for input qh, qk, and ql vectors for use with Horace
  
         Parameters
         ----------
         qh, qk, ql : (n_pts,) float ndarray
             The q-points to calculate at as separate vectors
-        pars : float ndarray
-            Parameters for the calculation (currently just the scale factor)
-            pars[0] = scale factor to multiply the intensity by
+        intensity_scale : float
+            The factor to multiply the intensity by
+        frequency_scale : float
+                The factor to multiply the phonon frequencies
+                by, as DFT can often overestimate frequencies
         args: tuple
             Arguments passed directly to the convolution function
         kwargs : dict
@@ -133,8 +134,6 @@ def horace_disp(self, qh, qk, ql, pars=[], *args, **kwargs):
         sf : (n_modes,) tuple of (n_pts,) float ndarray
             The dynamical structure corresponding to phonon energies in w as a tuple of numpy float vectors
         """
-
-        scale = (pars[0] if len(pars) > 0 else 1.) if hasattr(pars, '__len__') else pars
         if self.chunk > 0:
             lqh = len(qh)
             for i in range(int(np.ceil(lqh / self.chunk))):
@@ -157,8 +156,10 @@ def horace_disp(self, qh, qk, ql, pars=[], *args, **kwargs):
             if self.conversion_mat is not None:
                 qpts = np.matmul(qpts, self.conversion_mat)
             w, sf = self._calculate_sf(qpts)
-        if scale != 1.:
-            sf *= scale
+        if frequency_scale != 1.:
+            w *= frequency_scale
+        if intensity_scale != 1.:
+            sf *= intensity_scale
         sf = np.minimum(sf, self.lim)
         # Splits into different dispersion surfaces (python tuple == matlab cell)
         # But the data must be contiguous in memory so we need to do a real tranpose (.T just changes strides)

test/euphonic_horace_test.py

@@ -7,7 +7,6 @@
 import euphonic
 
 import euphonic_sqw_models
-
 def get_abspath(filename, sub_dir):
     test_dir = os.path.dirname(os.path.realpath(__file__))
     return test_dir + os.path.sep + sub_dir + os.path.sep + filename
@@ -23,9 +22,6 @@ def get_abspath(filename, sub_dir):
 conversion_mat = [None, (1./2)*np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]])]
 lim = [None, 1e2]  # Units: mbarn
 
-run_pars = [{'use_c': False, 'n_threads': 1, 'chunk': 5, 'dipole_parameter': 0.75},
-            {'use_c': True, 'n_threads': 1, 'chunk': 0},
-            {'use_c': True, 'n_threads': 2, 'chunk': 0}]
 
 qpts = [[0.0,  0.0,  0.0],
         [0.1,  0.2,  0.3],
@@ -36,29 +32,34 @@ def get_abspath(filename, sub_dir):
         [0.0, -0.5,  0.0],
         [0.0,  0.0, -0.5],
         [1.0, -1.0, -1.0]]
+qpts = np.array(qpts)
 scattering_lengths = {'La': 8.24, 'Zr': 7.16, 'O': 5.803,
                       'Si': 4.1491, 'Na': 3.63, 'Cl': 9.577}
-scale = 1.0
-qpts = np.array(qpts)
+iscale = 1.0
+freq_scale = 1.0
+pars = [iscale, freq_scale]
+
 
-def parameter_generator():
+def get_test_opts():
+    opts = []
     for tt in temp:
-        for mat in materials:
-            for dw in dw_grid:
-                for bos in bose:
-                    for ne in negative_e:
-                        for cmat in conversion_mat:
-                            for lm in lim:
-                                yield {'temperature': tt,
-                                       'material': mat,
-                                       'debye_waller_grid': dw,
-                                       'bose': bos,
-                                       'negative_e': ne,
-                                       'conversion_mat': cmat,
-                                       'lim': lm}
-
-def get_expected_output_filename(material_name, pars, opts):
-    fname = f"{material_name}_T{pars[0]}"
+        for dw in dw_grid:
+            for bos in bose:
+                for ne in negative_e:
+                    for cmat in conversion_mat:
+                        for lm in lim:
+                            opts.append({
+                                'temperature': tt,
+                                'debye_waller_grid': dw,
+                                'bose': bos,
+                                'negative_e': ne,
+                                'conversion_mat': cmat,
+                                'lim': lm})
+    return opts
+
+
+def get_expected_output_filename(material_name, opts):
+    fname = f"{material_name}_T{opts['temperature']}"
     if opts.get('debye_waller_grid', None) is not None:
         fname += "_dw" + "".join([str(v) for v in opts['debye_waller_grid']])
     if opts.get('bose', None) is not None:
@@ -80,6 +81,7 @@ def get_expected_output_filename(material_name, pars, opts):
     fname += '.mat'
     return get_abspath(fname, 'expected_output');
 
+
 def sum_degenerate_modes(w, sf):
     tol = 0.1;
     rows, cols = np.shape(w);
@@ -94,25 +96,29 @@ def sum_degenerate_modes(w, sf):
         summed_sf[i, :len(summed)] = summed;
     return summed_sf
 
-def calculate_w_sf(material_pars, material_constructor, par_dict):
-    fc = material_constructor(**material_pars)
-    par_dict['asr'] = 'reciprocal'
-    par_dict['scattering_lengths'] = scattering_lengths
-    coherent_sqw = euphonic_sqw_models.CoherentCrystal(fc, **par_dict)
-    w, sf = coherent_sqw.horace_disp(qpts[:,0], qpts[:,1], qpts[:,2], scale)
+
+def calculate_w_sf(material_opts, material_constructor, opt_dict,
+                   sum_sf=True, hdisp_args=(), hdisp_kwargs={}):
+    fc = material_constructor(**material_opts)
+    opt_dict['asr'] = 'reciprocal'
+    opt_dict['scattering_lengths'] = scattering_lengths
+    coherent_sqw = euphonic_sqw_models.CoherentCrystal(fc, **opt_dict)
+    w, sf = coherent_sqw.horace_disp(
+        qpts[:,0], qpts[:,1], qpts[:,2], *hdisp_args, **hdisp_kwargs)
     w = np.array(w).T
     sf = np.array(sf).T
     # Ignore acoustic structure factors by setting to zero - their
     # values can be unstable at small frequencies
     n = int(np.shape(sf)[1] / 2)
     sf[:,:3] = 0.
-    if 'negative_e' in par_dict and par_dict['negative_e']:
+    if 'negative_e' in opt_dict and opt_dict['negative_e']:
         sf[:,n:(n+3)] = 0
-    # Need to sum over degenerate modes to compare structure factors
-    sf_summed = sum_degenerate_modes(w, sf)
-    return w, sf_summed
+    if sum_sf:
+        sf = sum_degenerate_modes(w, sf)
+    return w, sf
 
-def get_expected_w_sf(fname):
+
+def get_expected_w_sf(fname, sum_sf=True):
     ref_dat = scipy.io.loadmat(fname)
     expected_w = np.concatenate(ref_dat['expected_w'][0], axis=1)
     expected_sf = np.concatenate(ref_dat['expected_sf'][0], axis=1)
@@ -121,51 +127,127 @@ def get_expected_w_sf(fname):
     if 'negative_e' in fname:
         n = int(np.shape(expected_sf)[1] / 2)
         expected_sf[:,n:(n+3)] = 0
+    if sum_sf:
+        expected_sf = sum_degenerate_modes(expected_w, expected_sf)
+    return expected_w, expected_sf
+
+@pytest.mark.parametrize("material", materials)
+@pytest.mark.parametrize("opt_dict", get_test_opts())
+@pytest.mark.parametrize("run_opts", [
+    {'use_c': False, 'n_threads': 1, 'chunk': 5, 'dipole_parameter': 0.75},
+    {'use_c': True, 'n_threads': 1, 'chunk': 0},
+    {'use_c': True, 'n_threads': 2}])
+def test_euphonic_sqw_models(material, opt_dict, run_opts):
+    material_name, material_constructor, material_opts = material
+    expected_w, expected_sf = get_expected_w_sf(
+        get_expected_output_filename(
+            material_name, opt_dict))
+
+    opt_dict.update(run_opts)
+    # Don't pass None options
+    opt_dict = {k: v for k, v in opt_dict.items() if v is not None}
+
+    w, sf = calculate_w_sf(material_opts, material_constructor, opt_dict)
+    npt.assert_allclose(w, expected_w, rtol=1e-5, atol=1e-2)
+    npt.assert_allclose(sf, expected_sf, rtol=1e-2, atol=1e-2)
+
+
+@pytest.mark.parametrize("material", materials)
+@pytest.mark.parametrize("opt_dict", [{
+    'temperature': 300,
+    'debye_waller_grid': [6, 6, 6],
+    'negative_e': True,
+    'conversion_mat': (1./2)*np.array([[-1, 1, 1],
+                                       [1, -1, 1],
+                                       [1, 1, -1]])}])
+# Test a combination of keyword and positional arguments
+# Is required for use with pace-python (Toby/Multifit
+# allow positional only for fitting parameters)
+@pytest.mark.parametrize("iscale, freqscale, args, kwargs",
+        [(1.0, 1.3, (), {'frequency_scale': 1.3}),
+         (1e-4, 1.0, (1e-4,), {}),
+         (1.5e3, 0.4, (1.5e3, 0.4), {}),
+         (2e3, 0.5, (), {'frequency_scale': 0.5, 'intensity_scale': 2e3}),
+         (50, 0.9, (50,), {'frequency_scale': 0.9})])
+def test_euphonic_sqw_models_pars(
+        material, opt_dict, iscale, freqscale, args, kwargs):
+    material_name, material_constructor, material_opts = material
+    expected_w, expected_sf = get_expected_w_sf(
+        get_expected_output_filename(
+            material_name, opt_dict), sum_sf=False)
+    w, sf = calculate_w_sf(material_opts, material_constructor, opt_dict,
+                           sum_sf=False, hdisp_args=args, hdisp_kwargs=kwargs)
+    # Sum sfs using the same frequencies - if expected_w and the scaled
+    # w are used, this could result in expected_sf and sf
+    # being summed differently
+    sf_summed = sum_degenerate_modes(expected_w, sf)
     expected_sf_summed = sum_degenerate_modes(expected_w, expected_sf)
-    return expected_w, expected_sf_summed
+    # Check that pars have scaled w and sf as expected
+    npt.assert_allclose(w, freqscale*expected_w, rtol=1e-5, atol=1e-2*freqscale)
+    npt.assert_allclose(sf_summed, iscale*expected_sf_summed,
+                        rtol=1e-2, atol=1e-2*iscale)
+
+
+@pytest.mark.parametrize("opt_dict", [{
+    'temperature': 300,
+    'debye_waller_grid': [6, 6, 6],
+    'negative_e': True,
+    'conversion_mat': (1./2)*np.array([[-1, 1, 1],
+                                       [1, -1, 1],
+                                       [1, 1, -1]])}])
+def test_old_behaviour_single_parameter_sets_intensity_scale(opt_dict):
+    iscale = 1.5
+    material_name, material_constructor, material_opts = materials[0]
+    fc = material_constructor(**material_opts)
+    opt_dict['asr'] = 'reciprocal'
+    opt_dict['scattering_lengths'] = scattering_lengths
+    coherent_sqw = euphonic_sqw_models.CoherentCrystal(fc, **opt_dict)
+    # Pass a single positional argument for iscale - this is the old
+    # syntax. Test that it hasn't broken, if it has we may need to raise
+    # a deprecation warning
+    w, sf = coherent_sqw.horace_disp(
+        qpts[:,0], qpts[:,1], qpts[:,2], [iscale])
+    w = np.array(w).T
+    sf = np.array(sf).T
+    # Ignore acoustic structure factors by setting to zero - their
+    # values can be unstable at small frequencies
+    n = int(np.shape(sf)[1] / 2)
+    sf[:,:3] = 0.
+    sf[:,n:(n+3)] = 0
 
-@pytest.mark.parametrize("par_dict", parameter_generator())
-def test_euphonic_sqw_models(par_dict):
-    material_name, material_constructor, material_pars = tuple(
-        par_dict.pop('material'))
     expected_w, expected_sf = get_expected_w_sf(
         get_expected_output_filename(
-            material_name, [par_dict['temperature'], scale], par_dict))
-
-    for run_par in run_pars:
-        par_dict.update(run_par)
-        for remove_if_none in list(par_dict.keys()):
-            if remove_if_none in par_dict and par_dict[remove_if_none] is None:
-                par_dict.pop(remove_if_none)
-        w, sf = calculate_w_sf(material_pars, material_constructor, par_dict)
-        npt.assert_allclose(w, expected_w, rtol=1e-5, atol=1e-2)
-        npt.assert_allclose(sf, expected_sf, rtol=1e-2, atol=1e-2)
-
-# For this test only test a subset of arguments
-def scale_test_parameter_generator():
-    for tt in temp:
-        for mat in materials:
-            yield {'temperature': tt,
-                   'material': mat,
-                   'debye_waller_grid':[6, 6, 6],
-                   'negative_e': True,
-                   'conversion_mat': (1./2)*np.array([[-1, 1, 1],
-                                                      [1, -1, 1],
-                                                      [1, 1, -1]])}
+            material_name, opt_dict), sum_sf=False)
+    # Sum sfs using the same frequencies - if expected_w and the scaled
+    # w are used, this could result in expected_sf and sf
+    # being summed differently
+    sf_summed = sum_degenerate_modes(expected_w, sf)
+    expected_sf_summed = sum_degenerate_modes(expected_w, expected_sf)
+    # Check that pars have scaled w and sf as expected
+    npt.assert_allclose(w, expected_w, rtol=1e-5, atol=1e-2)
+    npt.assert_allclose(sf_summed, iscale*expected_sf_summed,
+                        rtol=1e-2, atol=1e-2)
+
 
 # Units of sf have changed, test the calculated and old expected
 # values are the same apart from a scale factor
-@pytest.mark.parametrize("par_dict", scale_test_parameter_generator())
-def test_sf_unit_change(par_dict):
-    material_name, material_constructor, material_pars = tuple(
-        par_dict.pop('material'))
+@pytest.mark.parametrize("material", materials)
+@pytest.mark.parametrize("opt_dict", [{
+    'temperature': 300,
+    'debye_waller_grid': [6, 6, 6],
+    'negative_e': True,
+    'conversion_mat': (1./2)*np.array([[-1, 1, 1],
+                                       [1, -1, 1],
+                                       [1, 1, -1]])}])
+def test_sf_unit_change(material, opt_dict):
+    material_name, material_constructor, material_opts = material
     fname = get_expected_output_filename(
-            material_name, [par_dict['temperature'], scale], par_dict)
+            material_name, opt_dict)
     expected_w, expected_sf = get_expected_w_sf(os.path.join(
         os.path.dirname(fname),
         f'old_units_{os.path.basename(fname)}'))
 
-    w, sf = calculate_w_sf(material_pars, material_constructor, par_dict)
+    w, sf = calculate_w_sf(material_opts, material_constructor, opt_dict)
 
     npt.assert_allclose(w, expected_w, rtol=1e-5, atol=1e-2)