GEP fits

gistfile1.py

#!/usr/bin/env python

import os
from matplotlib import pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from exparser import Fitting

def fitGEPPoly(x, y, thr=90, deg=2, res=100, plot=True):

    """
    desc:
        Fit GEPs with a polynomial.

    arguments:
        x:  X data to fit
        y:  Y data to fit

    keywords:
        thr:    The threshold for the GEP
        deg:    The degrees for the polynomial, where 1 is a linear function.
        res:    The resolution of the fitted function.
        plot:   Whether the results should be plotted for debugging purposes.

    returns:
        A (x, y, xThr) tuple, where x and y are arrays with fitted data, and
        xThr is the first moment at which the y array exceeds the threshold\
        (thr).
    """

    p = np.polyfit(x, y, deg)
    fx = np.linspace(x.min(), x.max(), res)
    fy = p[-1]
    for i in range(1, deg+1):
        fy += p[-1-i]*fx**i
    i = np.where(fy > thr)
    if len(i) == 0:
        xThr = None
    else:
        xThr = fx[i[0][0]]
    if plot:
        plt.title('xThreshold: %s' % xThr)
        if xThr is not None:
            plt.axvline(xThr, linestyle=':', color='black')
        plt.axhline(thr, linestyle=':', color='black')
        plt.plot(x, y, 'o', color='black')
        plt.plot(fx, fy, ':', color='black')
        plt.show()
    return fx, fy, xThr

def fitGEPSigm(x, y, thr=90, res=100, plot=True):

    """
    desc:
        Fit GEPs with a sigmoid.

    arguments:
        x:  X data to fit
        y:  Y data to fit

    keywords:
        thr:    The threshold for the GEP
        res:    The resolution of the fitted function.
        plot:   Whether the results should be plotted for debugging purposes.

    returns:
        A (x, y, xThr) tuple, where x and y are arrays with fitted data, and
        xThr is the first moment at which the y array exceeds the threshold\
        (thr).
    """

    y = (y-50)/50
    try:
        opt, pcov = curve_fit(Fitting.sigmoid, x, y)
    except:
        opt = None
        xThr = None
        fy = None
        fx = None
    if opt is not None:
        fx = np.linspace(x.min(), x.max(), res)
        fy = Fitting.sigmoid(fx, x0=opt[0], k=opt[1])
        y = 50+50*y
        fy = 50+50*fy
        i = np.where(fy > thr)
        if len(i[0]) == 0:
            xThr = None
        else:
            xThr = fx[i[0][0]]
    if plot:
        plt.title('xThreshold: %s' % xThr)
        if xThr is not None:
            plt.axvline(xThr, linestyle=':', color='black')
        plt.axhline(thr, linestyle=':', color='black')
        plt.plot(x, y, 'o', color='black')
        if fy is not None:
            plt.plot(fx, fy, ':', color='black')
        plt.ylim(50, 100)
        plt.show()
    return fx, fy, xThr

for fname in sorted(os.listdir('npy')):
    if not fname.endswith('.npy'):
        continue
    path = 'npy/'+fname
    y = np.load(path)
    x = np.linspace(3, 10, 8)
    fx, fy, xThrPoly = fitGEPPoly(x, y, deg=2, plot=False)
    fx, fy, xThrSigm = fitGEPSigm(x, y, plot=False)
    plt.plot(xThrPoly, xThrSigm, 'o')
    print(fname, xThrPoly, xThrSigm)
plt.show()