Implement ODR fit

fitting.py

 import scipy
+import scipy.odr
 import numpy as np
 
 
@@ -31,12 +32,26 @@ def fit_lorenzian(S, V, log=False, **kwargs):
     return param, param_error, function, [X, function(X, *param)]
     
 
-def fit_function(function, x, y, p0=None, extend=0, **kwargs):
+def fit_function(function, x, y, *, xerr=None, yerr=None, p0=None, odr=None, extend=0, **kwargs):
     """
     :param log: if the y-data is log-scaled
     """
-    param, cov = scipy.optimize.curve_fit(function, x, y, p0, **kwargs)
-    param_error = np.sqrt(np.abs(cov.diagonal()))
+    if odr is None:
+        odr = xerr is not None or yerr is not None
+    
+    if odr:
+        # orthogonal distance regression
+        data = scipy.odr.RealData(x, y, sx=xerr, sy=yerr)
+        model = scipy.odr.Model(lambda beta, x: function(x, *beta))
+        odr = scipy.odr.ODR(data, model, p0)
+        # if not odr: odr.set_job(fit_type=2) # ordinary least-squares
+        output = odr.run()
+        param, param_error = output.beta, output.sd_beta
+    else:
+        # non-linear least squares
+        param, cov = scipy.optimize.curve_fit(function, x, y, p0, **kwargs)
+        param_error = np.sqrt(np.abs(cov.diagonal()))
+    
     mi, ma = min(x), max(x)
     X = np.linspace(mi - (ma-mi)*extend, ma + (ma-mi)*extend, 1000)
     return param, param_error, function, [X, function(X, *param)]