add regression citt, add new datafields to citt

src/paveit/postprocessing/citt.py

@@ -0,0 +1,131 @@
+import datetime
+
+import numpy as np
+from bson import ObjectId
+from paveit.datamodels import CITTSiffnessResults, RegCITT
+import lmfit as lm
+import pandas as pd
+import scipy.special as sf
+from scipy.optimize import curve_fit
+
+def temp_freq_equivalence(T, f, phi, T0=20.0):
+    
+    alphaT = np.exp(phi * ((1 / (T + 273.15)) - (1 / (T0 + 273.15))))
+    x = np.log(f * alphaT) / np.log(10)
+    
+    return x
+
+def stiffness_tp26(T, f, phi, Emax, Emin, z0, z1, T0=20.0):
+
+    x = temp_freq_equivalence(T, f, phi, T0)    
+    
+    E = Emin + (Emax - Emin) / (1 + np.exp(z1 * x + z0))
+
+    return E
+
+
+def calc_nu(T):
+    #TODO: Prüfen ob Formel stimmt!
+    nu = 0.15 + (0.35) / (1 + np.exp(3.1849 - 0.04233 * (9 / 5 * T + 32)))
+    return nu
+
+
+def citt(task_id: str):
+    """
+    Postprocessing task
+    """
+    
+    print('postprocessing')
+    
+    task_id = ObjectId(task_id)
+    
+    # read all data
+    data = []
+    parlist = ['f_set', 'T_set', 'stiffness', 'phase']
+    
+    for obj in CITTSiffnessResults.objects(task_id=task_id).only(*parlist):
+        data.append(dict((k, obj[k]) for k in parlist ))
+        
+    data = pd.DataFrame.from_records(data)
+    
+    #Emax/Emin
+    line_mod = lm.models.LinearModel()
+    out = line_mod.fit(data.stiffness, x=data.phase)
+    
+    print(out.best_values)
+    
+    Emax = line_mod.eval(out.params, x=0.0)
+    
+    Emin = 0
+    
+    assert Emin < Emax
+    
+    print(data.head())
+    
+    
+    # Fit data        
+    mod = lm.models.Model(stiffness_tp26, independent_vars=['f','T'])
+    
+    mod.set_param_hint(
+        'Emin',
+        value=Emin,
+        min=0,
+        max=0.9*Emax,
+        vary=True,
+    )
+    
+    mod.set_param_hint(
+        'Emax',
+        value=Emax,
+        min=0.9*Emax,
+        max=1.1*Emax,
+        vary=True,
+    )
+    
+    mod.set_param_hint(
+        'T0',
+        value=20.0,
+        vary=False,
+    )
+    
+    mod.set_param_hint('phi', value=25000, min=15000, max=35000, vary=True)
+    mod.set_param_hint('z0', value=1,min=1e-10, max=1000., vary=True)
+    mod.set_param_hint('z1', value=-1, min=-1000., max=-1e-10, vary=True)
+    
+    parms_fit = [
+        mod.param_hints['Emin']['value'], mod.param_hints['Emax']['value'],
+        mod.param_hints['phi']['value'], mod.param_hints['z0']['value'],
+        mod.param_hints['z1']['value']
+    ]
+    
+    
+    ## run fit
+    results = []
+    r2 = []
+
+    methods = ['leastsq', 'powell']
+
+    for method in methods:
+        result = mod.fit(data.stiffness, T=data.T_set, f=data.f_set, method=method, verbose=False)
+        
+        r2temp =  1.0 - result.redchi / np.var(data.stiffness.values, ddof=2)
+        r2.append(r2temp)
+
+        results.append(result)
+
+    best = np.nanargmax(r2)
+    
+    res = results[best].best_values
+    res['nsamples'] = len(data)
+    res['task_id'] = task_id
+    
+    res['stat_r2'] = r2[best]
+    
+    res['date'] = datetime.datetime.now()
+
+    print(res)
+    
+    # save results to db
+    doc = RegCITT.objects(task_id=task_id).modify(upsert=True, **res)
+    
+    return True