Python numpy.polyint函数代码示例

  def measureBdRatefct(self, reference, processed):
    """
    BJONTEGAARD    Bjontegaard metric calculation
    Bjontegaard's metric allows to compute the average % saving in bitrate
    between two rate-distortion curves [1].
    R1,Q1 - RD points for curve 1
    R2,Q2 - RD points for curve 2
    adapted from code from: (c) 2010 Giuseppe Valenzise
    """
    # numpy plays games with its exported functions.
    # pylint: disable=no-member
    # pylint: disable=too-many-locals
    # pylint: disable=bad-builtin
    R1 = [float(x[prX]) for x in reference]
    Q1 = [float(x[prY]) for x in reference]
    R2 = [float(x[prX]) for x in processed]
    Q2 = [float(x[prY]) for x in processed]

    #print(R1)
    #print(Q1)
    #print(R2)
    #print(Q2)

    log_R1 = map(math.log, R1)
    log_R2 = map(math.log, R2)

    log_R1 = numpy.log(R1)
    log_R2 = numpy.log(R2)

    #print(log_R1)
    #print(log_R2)

    # Best cubic poly fit for graph represented by log_ratex, psrn_x.
    poly1 = numpy.polyfit(Q1, log_R1, 3)
    poly2 = numpy.polyfit(Q2, log_R2, 3)

    # Integration interval.
    min_int = max([min(Q1), min(Q2)])
    max_int = min([max(Q1), max(Q2)])

    # find integral
    p_int1 = numpy.polyint(poly1)
    p_int2 = numpy.polyint(poly2)

    # Calculate the integrated value over the interval we care about.
    int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
    int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)

    # Calculate the average improvement.
    avg_exp_diff = (int2 - int1) / (max_int - min_int)

    # In really bad formed data the exponent can grow too large.
    # clamp it.
    if avg_exp_diff > 200:
      avg_exp_diff = 200

    # Convert to a percentage.
    avg_diff = (math.exp(avg_exp_diff) - 1) * 100

    return avg_diff

 def test_polyint_type(self) :
     """Ticket #944"""
     msg = "Wrong type, should be complex"
     x = np.ones(3, dtype=np.complex)
     assert_(np.polyint(x).dtype == np.complex, msg)
     msg = "Wrong type, should be float"
     x = np.ones(3, dtype=np.int)
     assert_(np.polyint(x).dtype == np.float, msg)

 def test_polyint_type(self):
     # Ticket #944
     msg = "Wrong type, should be complex"
     x = np.ones(3, dtype=complex)
     assert_(np.polyint(x).dtype == complex, msg)
     msg = "Wrong type, should be float"
     x = np.ones(3, dtype=int)
     assert_(np.polyint(x).dtype == float, msg)

def BdRate(group1, group2):
  """Compute the BD-rate between two score groups.

  The returned object also contains the range of PSNR values used
  to compute the result.

  Bjontegaard's metric allows to compute the average % saving in bitrate
  between two rate-distortion curves [1].

  rate1,psnr1 - RD points for curve 1
  rate2,psnr2 - RD points for curve 2

  adapted from code from: (c) 2010 Giuseppe Valenzise
  copied from code by [email protected], [email protected]

  """
  # pylint: disable=too-many-locals
  metric_set1 = group1.dataPoints()
  metric_set2 = group2.dataPoints()

  # numpy plays games with its exported functions.
  # pylint: disable=no-member
  # pylint: disable=bad-builtin
  psnr1 = [x[1] for x in metric_set1]
  psnr2 = [x[1] for x in metric_set2]

  log_rate1 = map(math.log, [x[0] for x in metric_set1])
  log_rate2 = map(math.log, [x[0] for x in metric_set2])

  # Best cubic poly fit for graph represented by log_ratex, psrn_x.
  poly1 = numpy.polyfit(psnr1, log_rate1, 3)
  poly2 = numpy.polyfit(psnr2, log_rate2, 3)

  # Integration interval.
  min_int = max([min(psnr1), min(psnr2)])
  max_int = min([max(psnr1), max(psnr2)])

  # find integral
  p_int1 = numpy.polyint(poly1)
  p_int2 = numpy.polyint(poly2)

  # Calculate the integrated value over the interval we care about.
  int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
  int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)

  # Calculate the average improvement.
  avg_exp_diff = (int2 - int1) / (max_int - min_int)

  # In really bad formed data the exponent can grow too large.
  # clamp it.
  if avg_exp_diff > 200:
    avg_exp_diff = 200

  # Convert to a percentage.
  avg_diff = (math.exp(avg_exp_diff) - 1) * 100

  return {'difference': avg_diff, 'psnr':[min_int, max_int]}

def bdrate(metric_set1, metric_set2):
  """
  BJONTEGAARD    Bjontegaard metric calculation
  Bjontegaard's metric allows to compute the average % saving in bitrate
  between two rate-distortion curves [1].

  rate1,psnr1 - RD points for curve 1
  rate2,psnr2 - RD points for curve 2

  adapted from code from: (c) 2010 Giuseppe Valenzise

  """
  rate1 = [x[0] for x in metric_set1]
  psnr1 = [x[1] for x in metric_set1]
  rate2 = [x[0] for x in metric_set2]
  psnr2 = [x[1] for x in metric_set2]

  log_rate1 = map(lambda x: math.log(x), rate1)
  log_rate2 = map(lambda x: math.log(x), rate2)

  # Best cubic poly fit for graph represented by log_ratex, psrn_x.
  p1 = numpy.polyfit(psnr1, log_rate1, 3)
  p2 = numpy.polyfit(psnr2, log_rate2, 3)

  # Integration interval.
  min_int = max([min(psnr1),min(psnr2)])
  max_int = min([max(psnr1),max(psnr2)])

  # find integral
  p_int1 = numpy.polyint(p1)
  p_int2 = numpy.polyint(p2)

  # Calculate the integrated value over the interval we care about.
  int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
  int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)

  # Calculate the average improvement.
  avg_exp_diff = (int2 - int1) / (max_int - min_int)

  # In really bad formed data the exponent can grow too large.
  # clamp it.
  if avg_exp_diff > 200 :
    avg_exp_diff = 200

  # Convert to a percentage.
  avg_diff = (math.exp(avg_exp_diff) - 1) * 100

  return avg_diff

def _sweep_poly_phase(t, poly):
    """
    Calculate the phase used by sweep_poly to generate its output.  See
    sweep_poly for a description of the arguments.

    """
    # polyint handles lists, ndarrays and instances of poly1d automatically.
    intpoly = polyint(poly)
    phase = 2*pi * polyval(intpoly, t)
    return phase

def BDPSNR(PSNR1, BR1, PSNR2, BR2):
    lBR1 = np.log10(BR1)
    p1 = np.polyfit( lBR1, PSNR1, 3)

    lBR2 = np.log10(BR2)
    p2 = np.polyfit( lBR2, PSNR2, 3)

    min_int = max(min(lBR1), min(lBR2))
    max_int = min(max(lBR1), max(lBR2))

    # find integral
    p_int1 = np.polyint(p1)
    p_int2 = np.polyint(p2)

    int1 = np.polyval(p_int1, max_int) - np.polyval(p_int1, min_int)
    int2 = np.polyval(p_int2, max_int) - np.polyval(p_int2, min_int)

    # find avg diff
    avg_diff = (int2-int1)/(max_int-min_int)

    return avg_diff

def bdsnr(metric_set1, metric_set2):
  """
  BJONTEGAARD    Bjontegaard metric calculation
  Bjontegaard's metric allows to compute the average gain in psnr between two
  rate-distortion curves [1].
  rate1,psnr1 - RD points for curve 1
  rate2,psnr2 - RD points for curve 2

  returns the calculated Bjontegaard metric 'dsnr'

  code adapted from code written by : (c) 2010 Giuseppe Valenzise
  http://www.mathworks.com/matlabcentral/fileexchange/27798-bjontegaard-metric/content/bjontegaard.m
  """
  rate1 = [x[0] for x in metric_set1]
  psnr1 = [x[1] for x in metric_set1]
  rate2 = [x[0] for x in metric_set2]
  psnr2 = [x[1] for x in metric_set2]

  log_rate1 = map(lambda x: math.log(x), rate1)
  log_rate2 = map(lambda x: math.log(x), rate2)

  # Best cubic poly fit for graph represented by log_ratex, psrn_x.
  p1 = numpy.polyfit(log_rate1, psnr1, 3)
  p2 = numpy.polyfit(log_rate2, psnr2, 3)

  # Integration interval.
  min_int = max([min(log_rate1),min(log_rate2)])
  max_int = min([max(log_rate1),max(log_rate2)])

  # Integrate p1, and p2.
  p_int1 = numpy.polyint(p1)
  p_int2 = numpy.polyint(p2)

  # Calculate the integrated value over the interval we care about.
  int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
  int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)

  # Calculate the average improvement.
  avg_diff = (int2 - int1) / (max_int - min_int)
  return avg_diff

def BDRate(PSNR1, BR1, PSNR2, BR2):
    lBR1 = np.log(BR1)
    p1 = np.polyfit( PSNR1, lBR1, 3)

    lBR2 = np.log(BR2)
    p2 = np.polyfit( PSNR2, lBR2, 3)

    min_int = max(min(PSNR1), min(PSNR2))
    max_int = min(max(PSNR1), max(PSNR2))

    # find integral
    p_int1 = np.polyint(p1)
    p_int2 = np.polyint(p2)

    int1 = np.polyval(p_int1, max_int) - np.polyval(p_int1, min_int)
    int2 = np.polyval(p_int2, max_int) - np.polyval(p_int2, min_int)

    # find avg diff
    avg_exp_diff = (int2-int1)/(max_int-min_int)
    avg_diff = (np.exp(avg_exp_diff)-1)*100

    return avg_diff

 def test_4(self):
   for type in classes:
     for M in range(type[1],type[2]+1):
       coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end)
       S = coll.Smat[1:,1:]
       # as in TEST 1, create and integrate a polynomial with random coefficients, but now of degree M-1
       poly_coeff = np.random.rand(M-1)
       poly_vals  = np.polyval(poly_coeff, coll.nodes)
       poly_int_coeff = np.polyint(poly_coeff)
       for i in range(1,M):
           int_ex = np.polyval(poly_int_coeff, coll.nodes[i]) - np.polyval(poly_int_coeff, coll.nodes[i-1])
           int_coll = np.dot(poly_vals, S[i,:])
           assert abs(int_ex - int_coll)<1e-12, "For node type " + type[0] + ", partial quadrature rule from Smat failed to integrate polynomial of degree M-1 exactly for M = " + str(M)

def chirp(t,f0=0,t1=1,f1=100,method='linear',phi=0,qshape=None):
    """Frequency-swept cosine generator.

    Inputs:

        t          --  array to evaluate waveform at
        f0, f1, t1 --  frequency (in Hz) of waveform is f0 at t=0 and f1 at t=t1
            Alternatively, if f0 is an array, then it forms the coefficients of
            a polynomial (c.f. numpy.polval()) in t. The values in f1, t1,
            method, and qshape are ignored.
        method     --  linear, quadratic, or logarithmic frequency sweep
        phi        --  optional phase in degrees
        qshape     --  shape parameter for quadratic curve: concave or convex
    """

    # Convert to radians.
    phi *= pi / 180
    if size(f0) > 1:
        # We were given a polynomial.
        return cos(2*pi*polyval(polyint(f0),t)+phi)
    if method in ['linear','lin','li']:
        beta = (f1-f0)/t1
        phase_angle = 2*pi * (f0*t + 0.5*beta*t*t)
    elif method in ['quadratic','quad','q']:
        if qshape == 'concave':
            mxf = max(f0,f1)
            mnf = min(f0,f1)
            f1,f0 = mxf, mnf
        elif qshape == 'convex':
            mxf = max(f0,f1)
            mnf = min(f0,f1)
            f1,f0 = mnf, mxf
        else:
            raise ValueError("qshape must be either 'concave' or 'convex' but "
                "a value of %r was given." % qshape)
        beta = (f1-f0)/t1/t1
        phase_angle = 2*pi * (f0*t + beta*t*t*t/3)
    elif method in ['logarithmic','log','lo']:
        if f1 <= f0:
            raise ValueError(
                "For a logarithmic sweep, f1=%f must be larger than f0=%f."
                % (f1, f0))
        beta = log10(f1-f0)/t1
        phase_angle = 2*pi * (f0*t + pow(10,beta*t)/(beta*log(10)))
    else:
        raise ValueError("method must be 'linear', 'quadratic', or "
            "'logarithmic' but a value of %r was given." % method)

    return cos(phase_angle + phi)

def f_evolution_element(x, y):
  root_real = 2.
  roots = np.zeros((3,3))
  if y < 0:
    dP = np.poly([root0, root_real + y * j, root_real - y * j])
  elif y > 0:
    dP = np.poly([root0, root_real+y, root_real-y])
  else:
    dP = np.poly([root0, root_real, -root_real])

  P = lamda*np.polyint(dP)
  cplx_roots = np.roots(dP)
  roots[:,0] = [_.real for _ in cplx_roots if _.real < max_x and _.real > min_x]
  roots[:,0] = np.sort(roots[:,0])
  z = np.polyval(P, x)
  for i in xrange(roots.shape[0]):
    roots[i,1] = y
    roots[i,2] = np.polyval(P, roots[i,0])
  return z,roots

def calc_omega(cp):
    cp.insert
    a=[]
    for i in range(len(cp)):
        ptmp = []
        tmp = 0
        for j in range(len(cp)):
            if j != i:
                row = []
                row.insert(0,1/(cp[i]-cp[j]))
                row.insert(1,-cp[j]/(cp[i]-cp[j]))
                ptmp.insert(tmp,row)
                tmp += 1
        p=[1]
        for j in range(len(cp)-1):
            p = conv(p,ptmp[j])
        pint = numpy.polyint(p)
        arow = []
        for j in range(len(cp)):
            arow.append(numpy.polyval(pint,cp[j]))
        a.append(arow)
    return a

  def test_1(self):
    for type in classes:
      for M in range(type[1],type[2]+1):
        coll = getattr(pySDC.CollocationClasses, type[0])(M, t_start, t_end)
        
        # some basic consistency tests
        assert np.size(coll.nodes)==np.size(coll.weights), "For node type " + type[0] + ", number of entries in nodes and weights is different"
        assert np.size(coll.nodes)==M, "For node type " + type[0] + ", requesting M nodes did not produce M entries in nodes and weights"


        # generate random set of polynomial coefficients
        poly_coeff = np.random.rand(coll.order-1)
        # evaluate polynomial at collocation nodes
        poly_vals  = np.polyval(poly_coeff, coll.nodes)
        # use python's polyint function to compute anti-derivative of polynomial
        poly_int_coeff = np.polyint(poly_coeff)
        # Compute integral from 0.0 to 1.0
        int_ex = np.polyval(poly_int_coeff, t_end) - np.polyval(poly_int_coeff, t_start)
        # use quadrature rule to compute integral
        int_coll = coll.evaluate(coll.weights, poly_vals)
        # For large values of M, substantial differences from different round of error have to be considered
        assert abs(int_ex - int_coll) < 1e-10, "For node type " + type[0] + ", failed to integrate polynomial of degree " + str(coll.order-1) + " exactly. Error: %5.3e" % abs(int_ex - int_coll)

def f_evolution(x, y):
  z = np.zeros((x.size, y.size))
  root_real = 2.
  roots = np.zeros((3,y.size,3))
  for k in xrange(y.size):
    if y[k] < 0:
      dP = np.poly([root0, root_real + y[k] * j, root_real - y[k] * j])
    elif y[k] > 0:
      dP = np.poly([root0, root_real + y[k], root_real-y[k]])
    else:
      dP = np.poly([root0, root_real, -root_real])

    P = lamda*np.polyint(dP)
    cplx_roots = np.roots(dP)
    roots[:,k,0] = [_.real for _ in cplx_roots if _.real < max_x and _.real > min_x]
    roots[:,k,0] = np.sort(roots[:,k,0])
    for i in xrange(x.size):
      z[i,k] = np.polyval(P, x[i])
    for i in xrange(roots.shape[0]):
      roots[i,k,1] = y[k]
      roots[i,k,2] = np.polyval(P, roots[i,k,0])
  return z,roots