init package and add first subpackages

src/paveit/analysis/__init__.py

@@ -0,0 +1,9 @@
+from .regression import *
+
+__all__ = [
+    # regession models
+    "fit_cos_simple",
+    "fit_cos",
+    #helper functions
+    "fit_cos_eval",
+]

src/paveit/analysis/regression.py

@@ -0,0 +1,159 @@
+import lmfit as lm
+import numpy as np
+import scipy.special as sf
+from scipy.optimize import curve_fit
+
+
+def cosfunc(t, A, w, p, c, e):
+    """
+    definition of the sin function
+    """
+
+    return A * np.cos(2 * np.pi * w * t + p) + e * t + c
+
+
+def fit_cos_eval(x, par):
+    """
+    
+    par: dict
+        fitting results
+    
+    """
+
+    ys = cosfunc(x, par['amp'], par['freq'], par['phase'], par['offset'],
+                 par['slope'])
+
+    return ys
+
+
+def regression_sine_fft():
+    """
+    fast fourier transformation for sine data
+    """
+
+    return []
+
+
+def fit_cos_simple(x, y, freq=10.0):
+    """
+    simple sine regression
+    
+    x: vector or list
+    y: vector or list
+    
+    freq: float
+        initial frequency of regression analysis
+        
+    
+    
+    """
+
+    tt = np.array(x)
+    yy = np.array(y)
+
+    guess_offset = np.mean(yy)
+    offset_b = 0.4 * abs(guess_offset)
+
+    guess_amp = abs(np.max(yy) - np.min(yy)) / 2.0
+
+    param_bounds = ([
+        0.3 * guess_amp, 0.99 * freq, -np.inf, guess_offset - offset_b, -np.inf
+    ], [1.3 * guess_amp, 1.01 * freq, np.inf, guess_offset + offset_b, np.inf])
+
+    popt, pcov = curve_fit(cosfunc, tt, yy, bounds=param_bounds)
+
+    A, w, p, c, e = popt
+
+    return {
+        "amp": A,
+        "freq": w,
+        "phase": p,
+        "offset": c,
+        "slope": e,
+    }
+
+
+def fit_cos(x, y, freq=10.0, constfreq=False):
+    """
+    sine regression
+    
+    x: vector or list
+    y: vector or list
+    
+    freq: float
+        initial frequency of regression analysis
+        
+    
+    
+    """
+
+    # step 1
+
+    res_step1 = fit_cos_simple(x, y, freq=freq)
+
+    # step 2: lmfit
+    mod = lm.models.Model(cosfunc)
+
+    mod.set_param_hint(
+        'A',
+        value=res_step1['amp'],
+        #min=res_step1['amp'] - 0.5 * abs(res_step1['amp']),
+        #max=res_step1['amp'] + 0.5 * abs(res_step1['amp'])
+    )
+
+    mod.set_param_hint(
+        'w',
+        value=freq,
+        vary=not constfreq,
+        #min=freq - 0.1 * freq,
+        #max=freq + 0.1 * freq,
+    )
+
+    mod.set_param_hint('p', value=res_step1['phase'], vary=True)
+
+    mod.set_param_hint('c', value=res_step1['offset'],
+                       vary=True)  #, min = -0.5, max = 0.5)
+
+    mod.set_param_hint('e', value=res_step1['slope'], vary=True)
+
+    parms_fit = [
+        mod.param_hints['A']['value'], mod.param_hints['w']['value'],
+        mod.param_hints['p']['value'], mod.param_hints['c']['value'],
+        mod.param_hints['e']['value']
+    ]
+
+    abweichung = []
+    chis = []
+    chis_red = []
+    results = []
+    r2 = []
+
+    methods = ['leastsq', 'powell']
+    dof = len(y) - len(parms_fit)
+
+    for method in methods:
+        result = mod.fit(y, t=x, method=method, verbose=False)
+        r2temp = 1 - result.residual.var() / np.var(y)
+        #    r2temp = result.redchi / np.var(yfit, ddof=2)
+        if r2temp < 0.:
+            r2temp = 0
+        r2.append(r2temp)
+
+        chi = result.chisqr
+        chis_red.append(result.redchi)
+        abweichung.append(sf.gammaincc(dof / 2., chi / 2))
+        chis.append(chi)
+        results.append(result)
+
+    res = {}
+    best = np.nanargmax(r2)
+
+    res[f'amp'] = results[best].best_values['A']
+    res[f'freq'] = results[best].best_values['w']
+    res[f'phase'] = results[best].best_values['p']
+    res[f'offset'] = results[best].best_values['c']
+    res[f'slope'] = results[best].best_values['e']
+
+    res[f'r2'] = r2[best]
+
+    return res