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

  def __init__(self):
    d = 3 # Degree of interpolating polynomial
    nk = 20 # Control discretization
    tf = 10.0 # End time
    h = tf/nk
    
    tau_root = [0] + collocation_points(d, "radau") # Choose collocation points
    C = NP.zeros((d+1,d+1)) # Coefficients of the collocation equation
    D = NP.zeros(d+1) # Coefficients of the continuity equation
    F = NP.zeros(d+1) # Coefficients of the quadrature function
    T = NP.zeros((nk,d+1)) # All collocation time points
    
    self.d = d
    self.nk = nk # Control discretization
    self.tf = tf # End time
    self.h = h # Size of the finite elements
    # Construct polynomial basis
    for j in range(d+1):
      # Construct Lagrange polynomials to get the polynomial basis at the collocation point
      p = NP.poly1d([1])
      for r in range(d+1):
        if r != j:
          p *= NP.poly1d([1, -tau_root[r]]) / (tau_root[j]-tau_root[r])
        D[j] = p(1.0) # Evaluate the polynomial at the final time to get the coefficients of the continuity equation
        pder = NP.polyder(p)
        for r in range(d+1): # Evaluate the time derivative of the polynomial at all collocation points to get the coefficients of the continuity equation
          C[j,r] = pder(tau_root[r])

          # Evaluate the integral of the polynomial to get the coefficients of the quadrature function
        pint = NP.polyint(p)
        F[j] = pint(1.0)
      for k in range(nk):
        for j in range(d+1):
          T[k,j] = h*(k + tau_root[j])
      
      self.T = T; self.C = C; self.D = D; self.F= F;

    def w(self, p_c, T = 293.15):
        """
        Moisture content [kg/m3]
        
        The capillary pressure p_c is required, and T is an optional argument
        """
        
        if self.w_method == 'vangenuchten':
        
            w = np.zeros(np.shape(p_c))
            n = 1./(1-self.w_m)
            for i in range(np.size(self.w_l)):
                w += self.w_sat * self.w_l[i] * \
                (1.+(self.w_alpha[i]*abs(p_c))**n[i])**(-self.w_m[i])
            
        elif self.w_method == 'polynomial':

            w = np.polyval(self.w_poly, ham.HR(p_c, T))
            
        elif self.w_method == 'slope':
            
            w = np.polyval(np.polyint(self.xi_poly), ham.HR(p_c, T))
            
        return w

if user_args.e:
    # Get experimental data
    directory = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
    filename = directory + "/experimental_results.txt"
    with open(filename) as input:
        data = zip(*(line.strip().split('\t') for line in input))
        data_name = data[0][0] + data[1][0] + data[2][0]
        exp_x = np.asfarray(data[0][1:])
        exp_y = np.asfarray(data[1][1:])
        # Error is given in % Rel Error
        exp_error = np.asfarray(data[2][1:])*exp_y/100.0

    print '\nLockwood Experimental'
    exp_fit = np.poly1d(np.polyfit(FMR, exp_y, poly))

    P = np.polyint(exp_fit)
    exp_int = P(FMR[-1])*CSDA_R
    print 'Integral using polyint      = ', exp_int

    x = np.linspace(0, FMR[-1], 1000)
    y = np.zeros(1000)
    for i in range(1000):
      y[i] = exp_fit(x[i])
    # plt.plot(x,y, color='b' )
    print 'Integral using simpson rule = ', integrate.simps(y, x)*CSDA_R

    # Plot the data
    line0,err0,arg3, = ax0.errorbar(FMR, exp_y, yerr=exp_error, label="Lockwood (Exp.)", fmt="-s", markersize=5 )


markers = ["--v","-.o",":^","--<","-.>",":+","--x","-.1",":2","--3","-.4",":8","--p","-.P",":*","--h","-.H",":X","--D","-.d"]

def getIntegral(p,interval):
    pint=np.polyint(p)
    return pint(interval[1])-pint(interval[0])

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

    Parameters
    ----------
    t : ndarray
        Times at which to evaluate the waveform.
    f0 : float or ndarray, optional
        Frequency (in Hz) of the waveform at time 0.  If `f0` is an
        ndarray, it specifies the frequency change as a polynomial in
        `t` (see Notes below).
    t1 : float, optional
        Time at which `f1` is specified.
    f1 : float, optional
        Frequency (in Hz) of the waveform at time `t1`.
    method : {'linear', 'quadratic', 'logarithmic'}, optional
        Kind of frequency sweep.
    phi : float
        Phase offset, in degrees.
    qshape : {'convex', 'concave'}
        If method is 'quadratic', `qshape` specifies its shape.

    Notes
    -----
    If `f0` is an array, it forms the coefficients of a polynomial in
    `t` (see `numpy.polval`). The polynomial determines the waveform
    frequency change in time.  In this case, the values of `f1`, `t1`,
    `method`, and `qshape` are ignored.

    This function is deprecated.  It will be removed in SciPy version 0.9.0.
    It exists so that during in version 0.8.0, the new chirp function can
    call this function to preserve the old behavior of the quadratic chirp.
    """
    warnings.warn("The function old_chirp is deprecated, and will be removed in "
                    "SciPy 0.9", DeprecationWarning)
    # 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)-1)/(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)