3
\$\begingroup\$

In my project I have to make curve-fitting with a lots of parameters, so scipy curve_fit struggles to find the answer.

For example: \$\ c_0 + c_1 \cdot cos (b_0 + b_1\cdot x + b_2\cdot x^2+ b_3\cdot x^3)\$ ,where \$ c_i, b_i \$ are the params to determine.

That's why I made a method which first tries to fit the desired function to only a little part of the data, then extends the area of fitting and uses the last optimal parameters as initial values for the next cycle. I wrote this object oriented, however I'm not sure it's nicely written or I don't know it's necessary to make it OO. I've never been taught OOP, this is what I've understood so far. The class has an attribute called obj. That's because it's embedded in a PyQt5 project, so I need the plot there. I modified the code to work on it's own. Here is an example dataset you can try.

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

########## these functions are actually imported from another file ###########
def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(value - array)).argmin()
    return array[idx], idx

def cos_fit1(x,c0, c1, b0, b1):
    return c0 + c1*np.cos(b0 + b1*x)

def cos_fit2(x,c0, c1, b0, b1, b2):
    return c0 + c1*np.cos(b0 + b1*x + b2*x**2)

def cos_fit3(x,c0, c1, b0, b1, b2, b3):
    return c0 + c1*np.cos(b0 + b1*x + b2*x**2 + b3*x**3)

def cos_fit4(x,c0, c1, b0, b1, b2, b3, b4):
    return c0 + c1*np.cos(b0 + b1*x + b2*x**2 + b3*x**3 + b4*x**4)

def cos_fit5(x,c0, c1, b0, b1, b2, b3, b4, b5):
    return c0 + c1*np.cos(b0 + b1*x + b2*x**2 + b3*x**3 + b4*x**4 + b5*x**5)
##############################################################################

class FitOptimizer(object):
    """Class to help achieve better fitting results."""
    def __init__(self, x, y, ref, sam, func=None, p0=None):
        self.x = x
        self.y = y
        self.ref = ref
        self.sam = sam
        if not isinstance(self.x, np.ndarray):
            try:
                self.x = np.asarray(self.x)
            except:
                raise
        if not isinstance(self.y, np.ndarray):
            try:
                self.y = np.asarray(self.y)
            except:
                raise
        if not isinstance(self.ref, np.ndarray):
            try:
                self.ref = np.asarray(self.ref)
            except:
                pass
        if not isinstance(self.sam, np.ndarray):
            try:
                self.sam = np.asarray(self.sam)
            except:
                pass
        if len(self.ref) == 0 or len(self.sam) == 0:
            self._y_norm = self.y
        else:
            self._y_norm = (self.y - self.ref - self.sam) / (2 *
                            np.sqrt(self.sam * self.ref))
        self.func = func
        if p0 is None:
            self.p0 = []
        else:
            self.p0 = p0
        self.popt = self.p0
        self._init_set = False
        self.obj = None
        self.counter = 0

    def set_initial_region(self, percent, center):
        """ Determines the initial region to fit"""
        self._init_set = True
        _, idx = find_nearest(self.x, center)
        self._upper_bound = np.floor(idx + (percent/2)*(len(self.x) + 1))
        self._lower_bound = np.floor(idx - (percent/2)*(len(self.x) + 1))
        self._upper_bound = self._upper_bound.astype(int)
        self._lower_bound = self._lower_bound.astype(int)
        if self._lower_bound < 0:
            self._lower_bound = 0
        if self._upper_bound > len(self.x):
            self._upper_bound = len(self.x)
        self._x_curr = self.x[self._lower_bound:self._upper_bound]
        self._y_curr = self._y_norm[self._lower_bound:self._upper_bound]

    def _extend_region(self, extend_by=0.2):
        """ Extends region of fit"""
        self._new_lower = np.floor(self._lower_bound - extend_by*len(self.x))
        self._new_upper = np.floor(self._upper_bound + extend_by*len(self.x))
        self._new_lower = self._new_lower.astype(int)
        self._new_upper = self._new_upper.astype(int)
        self._lower_bound = self._new_lower
        self._upper_bound = self._new_upper
        if self._new_lower < 0:
            self._new_lower = 0
        if self._new_upper > len(self.x):
            self._new_upper = len(self.x)
        self._x_curr = self.x[self._new_lower:self._new_upper]
        self._y_curr = self._y_norm[self._new_lower:self._new_upper]

    def _make_fit(self):
        """ Makes fit """
        try:
            if len(self._x_curr) == len(self.x):
                return True
            self.popt, self.pcov = curve_fit(self.func, self._x_curr, self._y_curr, 
                                             maxfev = 200000, p0 = self.p0)
            self.p0 = self.popt 
        except RuntimeError:
            if len(self.popt) > 4:
                self.p0[:3] = self.popt[:3] + np.random.normal(0, 100, len(self.popt)-3)
            else:
                self.p0 = self.popt + np.random.normal(0,100, len(self.popt))
            self.popt, self.pcov = curve_fit(self.func, self._x_curr, self._y_curr,
                                             maxfev = 200000, p0 = self.p0)

    def _fit_goodness(self):
        """ r^2 value"""
        residuals = self._y_curr - self.func(self._x_curr, *self.popt)
        ss_res = np.sum(residuals**2)
        ss_tot = np.sum((self._y_curr - np.mean(self._y_curr))**2)
        return 1 - (ss_res / ss_tot)

    def _show_fit(self):
        """ Shows fit on self.obj """
        try:
            self.obj.plot(self._x_curr, self._y_curr, 'k-', label = 'Affected data')
            self.obj.plot(self._x_curr, self.func(self._x_curr, *self.popt),
                               'r--', label = 'Fit')
            self.obj.show()
        except Exception as e:
            print(e)


    def run(self, r_extend_by, r_threshold=0.9, max_tries=10000):
        if self._init_set == False:
            raise ValueError('Set the initial conditions.')
        self._make_fit()
        while self._fit_goodness() > r_threshold:
            self._extend_region(r_extend_by)
            self._make_fit()
            self.counter +=1
            if self._make_fit() == True:
                self._show_fit()
                return self.popt
                break
            if self.counter == max_tries:
                self._show_fit()
                return np.zeros_like(self.popt)
                break

        while self._fit_goodness() < r_threshold:
            self._make_fit()
            self.counter +=1
            if self.counter == max_tries:
                self._show_fit()
                return np.zeros_like(self.popt)
                break


""" EXAMPLE """
a, b, c, d = np.loadtxt('test.txt', delimiter = ',', unpack = True )
f = FitOptimizer(a, b, c, d, func = cos_fit3)
f.obj = plt
f.p0 = [1,1,1,1,1,1]
f.set_initial_region(0.2, 2.4)
f.run(r_extend_by = 0.1, r_threshold = 0.85)

I appreciate any improvements in the code.

\$\endgroup\$
1
\$\begingroup\$

Don't repeat yourself

Your family of cos_fit functions can call the full-form function, e.g.

def cos_fit1(x,c0, c1, b0, b1):
    return cos_fit5(x, c0, c1, b0, b1, 0, 0, 0, 0)

You may also want to consider using built-in polynomial support, e.g. https://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.polyval.html#numpy.polynomial.polynomial.polyval

As for these expressions:

    self._upper_bound = np.floor(idx + (percent/2)*(len(self.x) + 1))
    self._lower_bound = np.floor(idx - (percent/2)*(len(self.x) + 1))

reuse the inner term:

delta = percent * (1 + len(self.x)) / 2
self._upper_bound = np.floor(idx + delta)
self._lower_bound = np.floor(idx - delta)

No-op except

Delete this:

        except:
            raise

It doesn't do anything helpful.

As for your exception-swallowing

        except:
            pass

this is a deeply bad idea. You need to narrow the exception type caught to the one you're actually expecting. If the program spends any length of time in the try block, this form of except will prevent Ctrl+C break from working.

Order of operations

(self.y - self.ref - self.sam) / (2 *
                            np.sqrt(self.sam * self.ref))

can be

(self.y - self.ref - self.sam) / 2 / np.sqrt(self.sam * self.ref)

Similarly, 1 - (ss_res / ss_tot) doesn't need parens.

This is a comment

def _make_fit(self):
    """ Makes fit """

No kidding! Either your comments should add more than the method name, or you should just delete them.

No-op break

            return np.zeros_like(self.popt)
            break

The break is never executed due to the return, so you can delete it.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.