Merge pull request #99 from igresh/delta_residual_fix

Delta residual fix

Author: andyfaff

refellips/dataSE.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-""""
+"""
 A basic representation of a 1D dataset
 """
 

refellips/objectiveSE.py

@@ -18,6 +18,28 @@
 from refnx._lib import flatten
 
 
+def circular_distance(angle1, angle2, period=360):
+    """
+    Calculates the circular distance between two angles.
+    Units (rad or deg) do not matter as long as they are
+    consistent.
+
+    Parameters
+    ----------
+    angle1 : float, np.ndarray
+        First angle
+    angle2 : float, np.ndarray
+        Second angle
+    period : float
+        The period of the circular domain (e.g.,360 for full circle).
+
+    Returns:
+        float or np.ndarray: The shortest circular distance between the angles.
+    """
+    diff = np.abs(angle1 - angle2)
+    return np.minimum(diff, period - diff)
+
+
 class ObjectiveSE(BaseObjective):
     """
     Objective function for using with curvefitters such as
@@ -178,9 +200,9 @@ def residuals(self, pvals=None):
 
         wavelength, aoi, psi_d, delta_d = self.data.data
         wavelength_aoi = np.c_[wavelength, aoi]
-
         psi, delta = self.model(wavelength_aoi)
-        return np.r_[psi - psi_d, delta - delta_d]
+        delta_err = circular_distance(delta, delta_d, period=360)
+        return np.r_[psi - psi_d, delta_err]
 
     def chisqr(self, pvals=None):
         """
@@ -369,16 +391,16 @@ def logl(self, pvals=None):
 
         psi, delta = self.model(wavelength_aoi)
 
-        model = np.r_[psi, delta]
-
         logl = 0.0
 
         # TODO investigate ellipsometry uncertainties
         # here just set it to unity
         y_err = 1
         if self.lnsigma is not None:
+            _model = np.r_[psi, delta]
             var_y = (
-                y_err * y_err + np.exp(2 * float(self.lnsigma)) * model * model
+                y_err * y_err
+                + np.exp(2 * float(self.lnsigma)) * _model * _model
             )
         else:
             var_y = y_err**2
@@ -387,7 +409,8 @@ def logl(self, pvals=None):
         if self.weighted:
             logl += np.log(2 * np.pi * var_y)
 
-        logl += (np.r_[psi_d, delta_d] - model) ** 2 / var_y
+        res = self.residuals(None)
+        logl += (res) ** 2 / var_y
 
         # nans play havoc
         if np.isnan(logl).any():

refellips/tests/test_refellips.py

@@ -477,3 +477,42 @@ def test_TaucLorentz():
 
     assert_allclose(psi, data.psi, rtol=0.0176)
     assert_allclose(delta, data.delta, rtol=0.0066)
+
+
+def test_residuals():
+    """
+    Specifically testing for correct delta error handling (delta is an angle between 0 and 360)
+    """
+    wavelength = np.linspace(400, 900, 100)
+    aoi = 70
+
+    si = load_material("silicon")
+    sio2 = load_material("silica")
+    polymer = Cauchy(A=1.43, B=0.0005)
+    solvent = Cauchy(A=1.45, B=0.0005)
+
+    polymer_layer_1 = polymer(55)
+    polymer_layer_2 = polymer(45)
+
+    struc_1 = solvent() | polymer_layer_1 | sio2(20) | si()
+    struc_2 = solvent() | polymer_layer_2 | sio2(20) | si()
+
+    model_1 = ReflectModelSE(struc_1)
+    model_2 = ReflectModelSE(struc_2)
+
+    psi_1, delta_1 = model_1(np.c_[wavelength, np.ones_like(wavelength) * aoi])
+    psi_2, delta_2 = model_2(np.c_[wavelength, np.ones_like(wavelength) * aoi])
+    faux_data = DataSE(
+        [wavelength, np.ones_like(wavelength) * aoi, psi_1, delta_1]
+    )
+
+    obj = ObjectiveSE(model=model_2, data=faux_data)
+    res = obj.residuals()
+    obj_psi_res = res[: int(len(res) / 2)]
+    obj_delta_res = res[int(len(res) / 2) :]
+
+    test_psi_res = psi_2 - psi_1
+    test_delta_res = np.abs((delta_2 - delta_1 + 180) % 360 - 180)
+
+    assert_allclose(test_psi_res, obj_psi_res, rtol=0.002)
+    assert_allclose(test_delta_res, obj_delta_res, rtol=0.003)