(points, x_min_bound=None, x_max_bound=None, plot=False)
return x_sol
Gives the minimizer and minimum of the interpolating polynomial over given points
    based on function and derivative information. Defaults to bisection if no critical
    points are valid.

    Based on polyinterp.m Matlab function in minFunc by Mark Schmidt with some slight
    modifications.

    Implemented by: Hao-Jun Michael Shi and Dheevatsa Mudigere
    Last edited 12/6/18.

    Inputs:
        points (nparray): two-dimensional array with each point of form [x f g]
        x_min_bound (float): minimum value that brackets minimum (default: minimum of points)
        x_max_bound (float): maximum value that brackets minimum (default: maximum of points)
        plot (bool): plot interpolating polynomial

    Outputs:
        x_sol (float): minimizer of interpolating polynomial
        F_min (float): minimum of interpolating polynomial

    Note:
      . Set f or g to np.nan if they are unknown
Gives the minimizer and minimum of the interpolating polynomial over given points
    based on function and derivative information. Defaults to bisection if no critical
    points are valid.
Gives the minimizer and minimum of the interpolating polynomial over given points based on function and derivative information . Defaults to bisection if no critical points are valid .
def polyinterp(points, x_min_bound=None, x_max_bound=None, plot=False):
    """
    Gives the minimizer and minimum of the interpolating polynomial over given points
    based on function and derivative information. Defaults to bisection if no critical
    points are valid.

    Based on polyinterp.m Matlab function in minFunc by Mark Schmidt with some slight
    modifications.

    Implemented by: Hao-Jun Michael Shi and Dheevatsa Mudigere
    Last edited 12/6/18.

    Inputs:
        points (nparray): two-dimensional array with each point of form [x f g]
        x_min_bound (float): minimum value that brackets minimum (default: minimum of points)
        x_max_bound (float): maximum value that brackets minimum (default: maximum of points)
        plot (bool): plot interpolating polynomial

    Outputs:
        x_sol (float): minimizer of interpolating polynomial
        F_min (float): minimum of interpolating polynomial

    Note:
      . Set f or g to np.nan if they are unknown

    """
    no_points = points.shape[0]
    order = np.sum(1 - np.isnan(points[:,1:3]).astype('int')) - 1

    x_min = np.min(points[:, 0])
    x_max = np.max(points[:, 0])

    # compute bounds of interpolation area
    if(x_min_bound is None):
        x_min_bound = x_min
    if(x_max_bound is None):
        x_max_bound = x_max

    # explicit formula for quadratic interpolation
    if no_points == 2 and order == 2 and plot is False:
        # Solution to quadratic interpolation is given by:
        # a = -(f1 - f2 - g1(x1 - x2))/(x1 - x2)^2
        # x_min = x1 - g1/(2a)
        # if x1 = 0, then is given by:
        # x_min = - (g1*x2^2)/(2(f2 - f1 - g1*x2))

        if(points[0, 0] == 0):
            x_sol = -points[0, 2]*points[1, 0]**2/(2*(points[1, 1] - points[0, 1] - points[0, 2]*points[1, 0]))
        else:
            a = -(points[0, 1] - points[1, 1] - points[0, 2]*(points[0, 0] - points[1, 0]))/(points[0, 0] - points[1, 0])**2
            x_sol = points[0, 0] - points[0, 2]/(2*a)

        x_sol = np.minimum(np.maximum(x_min_bound, x_sol), x_max_bound)

    # explicit formula for cubic interpolation
    elif no_points == 2 and order == 3 and plot is False:
        # Solution to cubic interpolation is given by:
        # d1 = g1 + g2 - 3((f1 - f2)/(x1 - x2))
        # d2 = sqrt(d1^2 - g1*g2)
        # x_min = x2 - (x2 - x1)*((g2 + d2 - d1)/(g2 - g1 + 2*d2))
        d1 = points[0, 2] + points[1, 2] - 3*((points[0, 1] - points[1, 1])/(points[0, 0] - points[1, 0]))
        d2 = np.sqrt(d1**2 - points[0, 2]*points[1, 2])
        if np.isreal(d2):
            x_sol = points[1, 0] - (points[1, 0] - points[0, 0])*((points[1, 2] + d2 - d1)/(points[1, 2] - points[0, 2] + 2*d2))
            x_sol = np.minimum(np.maximum(x_min_bound, x_sol), x_max_bound)
        else:
            x_sol = (x_max_bound + x_min_bound)/2

    # solve linear system
    else:
        # define linear constraints
        A = np.zeros((0, order+1))
        b = np.zeros((0, 1))

        # add linear constraints on function values
        for i in range(no_points):
            if not np.isnan(points[i, 1]):
                constraint = np.zeros((1, order+1))
                for j in range(order, -1, -1):
                    constraint[0, order - j] = points[i, 0]**j
                A = np.append(A, constraint, 0)
                b = np.append(b, points[i, 1])

        # add linear constraints on gradient values
        for i in range(no_points):
            if not np.isnan(points[i, 2]):
                constraint = np.zeros((1, order+1))
                for j in range(order):
                    constraint[0, j] = (order-j)*points[i,0]**(order-j-1)
                A = np.append(A, constraint, 0)
                b = np.append(b, points[i, 2])

        # check if system is solvable
        if(A.shape[0] != A.shape[1] or np.linalg.matrix_rank(A) != A.shape[0]):
            x_sol = (x_min_bound + x_max_bound)/2
            f_min = np.Inf
        else:
            # solve linear system for interpolating polynomial
            coeff = np.linalg.solve(A, b)

            # compute critical points
            dcoeff = np.zeros(order)
            for i in range(len(coeff) - 1):
                dcoeff[i] = coeff[i]*(order-i)

            crit_pts = np.array([x_min_bound, x_max_bound])
            crit_pts = np.append(crit_pts, points[:, 0])

            if not np.isinf(dcoeff).any():
                roots = np.roots(dcoeff)
                crit_pts = np.append(crit_pts, roots)

            # test critical points
            f_min = np.Inf
            x_sol = (x_min_bound + x_max_bound)/2 # defaults to bisection
            for crit_pt in crit_pts:
                if np.isreal(crit_pt) and crit_pt >= x_min_bound and crit_pt <= x_max_bound:
                    F_cp = np.polyval(coeff, crit_pt)
                    if np.isreal(F_cp) and F_cp < f_min:
                        x_sol = np.real(crit_pt)
                        f_min = np.real(F_cp)

            if (plot):
                import matplotlib.pyplot as plt
                plt.figure()
                x = np.arange(x_min_bound, x_max_bound, (x_max_bound - x_min_bound)/10000)
                f = np.polyval(coeff, x)
                plt.plot(x, f)
                plt.plot(x_sol, f_min, 'x')

    return x_sol
def polyinterp ( points , x_min_bound = None , x_max_bound = None , plot = False ) : no_points = points . shape [ 0 ] order = np . sum ( 1 - np . isnan ( points [ : , 1 : 3 ] ) . astype ( 'int' ) ) - 1 x_min = np . min ( points [ : , 0 ] ) x_max = np . max ( points [ : , 0 ] ) # compute bounds of interpolation area if ( x_min_bound is None ) : x_min_bound = x_min if ( x_max_bound is None ) : x_max_bound = x_max # explicit formula for quadratic interpolation if no_points == 2 and order == 2 and plot is False : # Solution to quadratic interpolation is given by: # a = -(f1 - f2 - g1(x1 - x2))/(x1 - x2)^2 # x_min = x1 - g1/(2a) # if x1 = 0, then is given by: # x_min = - (g1*x2^2)/(2(f2 - f1 - g1*x2)) if ( points [ 0 , 0 ] == 0 ) : x_sol = - points [ 0 , 2 ] * points [ 1 , 0 ] ** 2 / ( 2 * ( points [ 1 , 1 ] - points [ 0 , 1 ] - points [ 0 , 2 ] * points [ 1 , 0 ] ) ) else : a = - ( points [ 0 , 1 ] - points [ 1 , 1 ] - points [ 0 , 2 ] * ( points [ 0 , 0 ] - points [ 1 , 0 ] ) ) / ( points [ 0 , 0 ] - points [ 1 , 0 ] ) ** 2 x_sol = points [ 0 , 0 ] - points [ 0 , 2 ] / ( 2 * a ) x_sol = np . minimum ( np . maximum ( x_min_bound , x_sol ) , x_max_bound ) # explicit formula for cubic interpolation elif no_points == 2 and order == 3 and plot is False : # Solution to cubic interpolation is given by: # d1 = g1 + g2 - 3((f1 - f2)/(x1 - x2)) # d2 = sqrt(d1^2 - g1*g2) # x_min = x2 - (x2 - x1)*((g2 + d2 - d1)/(g2 - g1 + 2*d2)) d1 = points [ 0 , 2 ] + points [ 1 , 2 ] - 3 * ( ( points [ 0 , 1 ] - points [ 1 , 1 ] ) / ( points [ 0 , 0 ] - points [ 1 , 0 ] ) ) d2 = np . sqrt ( d1 ** 2 - points [ 0 , 2 ] * points [ 1 , 2 ] ) if np . isreal ( d2 ) : x_sol = points [ 1 , 0 ] - ( points [ 1 , 0 ] - points [ 0 , 0 ] ) * ( ( points [ 1 , 2 ] + d2 - d1 ) / ( points [ 1 , 2 ] - points [ 0 , 2 ] + 2 * d2 ) ) x_sol = np . minimum ( np . maximum ( x_min_bound , x_sol ) , x_max_bound ) else : x_sol = ( x_max_bound + x_min_bound ) / 2 # solve linear system else : # define linear constraints A = np . zeros ( ( 0 , order + 1 ) ) b = np . zeros ( ( 0 , 1 ) ) # add linear constraints on function values for i in range ( no_points ) : if not np . isnan ( points [ i , 1 ] ) : constraint = np . zeros ( ( 1 , order + 1 ) ) for j in range ( order , - 1 , - 1 ) : constraint [ 0 , order - j ] = points [ i , 0 ] ** j A = np . append ( A , constraint , 0 ) b = np . append ( b , points [ i , 1 ] ) # add linear constraints on gradient values for i in range ( no_points ) : if not np . isnan ( points [ i , 2 ] ) : constraint = np . zeros ( ( 1 , order + 1 ) ) for j in range ( order ) : constraint [ 0 , j ] = ( order - j ) * points [ i , 0 ] ** ( order - j - 1 ) A = np . append ( A , constraint , 0 ) b = np . append ( b , points [ i , 2 ] ) # check if system is solvable if ( A . shape [ 0 ] != A . shape [ 1 ] or np . linalg . matrix_rank ( A ) != A . shape [ 0 ] ) : x_sol = ( x_min_bound + x_max_bound ) / 2 f_min = np . Inf else : # solve linear system for interpolating polynomial coeff = np . linalg . solve ( A , b ) # compute critical points dcoeff = np . zeros ( order ) for i in range ( len ( coeff ) - 1 ) : dcoeff [ i ] = coeff [ i ] * ( order - i ) crit_pts = np . array ( [ x_min_bound , x_max_bound ] ) crit_pts = np . append ( crit_pts , points [ : , 0 ] ) if not np . isinf ( dcoeff ) . any ( ) : roots = np . roots ( dcoeff ) crit_pts = np . append ( crit_pts , roots ) # test critical points f_min = np . Inf x_sol = ( x_min_bound + x_max_bound ) / 2 # defaults to bisection for crit_pt in crit_pts : if np . isreal ( crit_pt ) and crit_pt >= x_min_bound and crit_pt <= x_max_bound : F_cp = np . polyval ( coeff , crit_pt ) if np . isreal ( F_cp ) and F_cp < f_min : x_sol = np . real ( crit_pt ) f_min = np . real ( F_cp ) if ( plot ) : import matplotlib . pyplot as plt plt . figure ( ) x = np . arange ( x_min_bound , x_max_bound , ( x_max_bound - x_min_bound ) / 10000 ) f = np . polyval ( coeff , x ) plt . plot ( x , f ) plt . plot ( x_sol , f_min , 'x' ) return x_sol