title: Python optimizing a calculation with scipy.integrate.quad (takes very long)
body: <p>I´m currently writing a program in python for calculating the total spectral emissivity (infrared waves) of any given material at different temperatures (200K - 500K), based on measurement data received by measuring the directional - hemispherical emissivity of the material at many different wavelengths using an IR - spectroscope. The calculation is done by integrating the measured intensity over all wavelenghts, using Plancks law as a weighing function (all of this doesn´t really matter for my question itself; i just want to explain the background so that the code is easier to understand). This is my code:</p>
  
  <pre><code>from scipy import integrate
  from scipy.interpolate import CubicSpline
  import numpy as np
  import math as m
  
  
  def planck_blackbody(lambda_, T):  # wavelength, temperature
      h = float(6.6260755e-34)
      c = float(2.99792458e+8)
      k = float(1.380658e-23)
      try:
          a = 2.0 * h * (c ** 2)
          b = h * c / (lambda_ * k * T)
          intensity = a / ((lambda_ ** 5) * (m.exp(b) - 1.0))
          return float(intensity)
      except OverflowError:  # for lower temperatures
          pass
  
  
  def spectral_emissivity(emifilename, t, lambda_1, lambda_2):
      results = []
  
      with open(emifilename, 'r') as emifile:
         emilines = emifile.readlines()
         try:
             w = [float(x.split('\t')[0].strip('\n')) * 1e-6 for x in emilines]
             e = [float(x.split('\t')[1].strip('\n')) for x in emilines]
         except ValueError:
             pass
      w = np.asarray(w)  # wavelength
      e = np.asarray(e)  # measured intensity
  
      def part_1(lambda_, T):
          E = interp1d(w, e, fill_value = 'extrapolate')(lambda_)
          return E * planck_blackbody(lambda_, T)
  
      def E_complete(T):
          E_complete_part_1 = integrate.quad(part_1, lambda_1, lambda_2, args=T, limit=50)
          E_complete_part_2 = integrate.quad(planck_blackbody, lambda_1, lambda_2, args=T, limit=50)
          return E_complete_part_1[0] / E_complete_part_2[0]
  
      for T in t:
          results.append([T, E_complete(T)])
  
      with open("{}.plk".format(emifilename[:-4]), 'w') as resultfile:
          for item in results:
              resultfile.write("{}\t{}\n".format(item[0], item[1]))
  
  t = np.arange(200, 501, 1)
  spectral_emissivity('C:\test.dat', t, 1.4e-6, 35e-6)
  </code></pre>
  
  <p>The measured intensity is stored in a text file with two columns, the first being the wavelength of the infrared waves and the second being the directional-hemispherical emissivity of the measured material at that wavelength.</p>
  
  <p>When i run this code, while it is producing the right results, i still encounter 2 problems:</p>
  
  <ol>
  <li><p>I get an error message from scipy.integrate.quad:</p>
  
  <pre><code>  IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
    If increasing the limit yields no improvement it is advised to analyze 
    the integrand in order to determine the difficulties.  If the position of a 
    local difficulty can be determined (singularity, discontinuity) one will 
    probably gain from splitting up the interval and calling the integrator 
    on the subranges.  Perhaps a special-purpose integrator should be used.
    warnings.warn(msg, IntegrationWarning)
  </code></pre></li>
  </ol>
  
  <p>Can someone explain to me what exactly this means? I understand that integrate.quad is a numerical iteration method and that my functions somehow seem to require more than 50 iterations, but is there a way around this? i tried increasing the Limit, but even with 200 i still get this error message... it´s especially weird given that the integrands are pretty straightforward functions...</p>
  
  <ol start="2">
  <li>is closely connected to the first problem: this program takes ages (about 5 minutes!) to finish one single file, but i need to process many files every hour. cProfile reveals that 98% percent of this time is spent inside the integraion function. A MathCad program doing the exact same thing and producing the same outputs only takes some seconds to finish. Even though i spent the last week seatching for a solution, i simply don´t manage to speed this program up, and no one else on stackoverflow and elsewhere seems to have comparable time problems with integrate.quad.</li>
  </ol>
  
  <p>So, finally, my question: is there any obvious way to optimize this code for it to run faster (except from compiling it into C+ or anything like that)? I tried reducing all floats to 6 digits (i can´t go any lower in accuracy) but that didn´t change anything.</p>
  
  <p><strong>Update:</strong> looking into it some more, i figured out that most of the time wasn´t actually consumed by the Integration itself, but by the <code>CubicSpline</code> operation that i used to interpolate my data. I tried out different methods and <code>CubicSpline</code> seemed to be the only working one for some reason (even though my data is monotonically increasing, i got errors from every other method i tried, saying that some values were either above or below the interpolation range). That is, until i found out about extrapolation <code>with scipy.interpolate.interp1d</code>and <code>(fill_value = 'extrapolate'</code>. Ths did the trick for me, enabling me to use the far less consuming <code>interp1d</code> method and effectively reducing the runtime of my program from 280 to 49 seconds (also added list comprehension for <code>w</code> and <code>e</code>). While this is a big improvement, i still wonder why my program takes nearly 1 Minute to calculate some integrals... and i still get the above mentioned <code>IntegrationWarning</code>. So any advice is highly appreciated!</p>
  
  <p>(btw, since i am pretty new to python, I´m happy about any tips or critique i can get!)</p>