# [Tutor] Fitting data to error function

Colin Ross colin.ross.dal at gmail.com
Sun Mar 15 16:27:50 CET 2015

Hi Danny,

Thanks for the help! As you mentioned, using scipy.special.erfc was a much
better idea. Below is a copy of my program and the stack trace, showing a
new error. It seems that the first auto correlation works,  however the
second fails.

###############################

# Autocorrelation program

###############################

import numpy as np

from numpy import ma, logical_or

import pylab

from pylab import *

import matplotlib.pyplot as plt

from matplotlib.ticker import AutoMinorLocator

import math

from math import exp

import scipy

import sys

from scipy.optimize import curve_fit

from scipy.fftpack import fft, ifft, fftfreq

##############################

#fitting function (gaussian)

##############################

def fit(x,a,b,c):

return a*np.exp(-(x-b)**2/c)

#(1.0/(np.sqrt(2*np.pi)*sigma))*np.exp(-(x-mu)**2 / (2*sigma**2))

##############################

# Load data from .txt file

##############################

data1 = np.loadtxt('TL_pulseC.txt',skiprows = 2 ,usecols = (0,1))

data2 = np.loadtxt('cos_pulseC.txt',skiprows = 2 ,usecols = (0,1))

#print "\n  Stage Position (um)     Amplitude (V)"

#print data

################################

# Create data arrays

################################

pos = np.array(data1[:,0])           # micrometers

pos_m = pos*1.0*10**(-6)             # meters

print pos_m

amp = np.array(data1[:,1])           # V

amp_auto = np.array(data2[:,1])

#################################################################################################

# Finding the index of the maximum amplitude where the pulse delay is zero.
Then calculating the

# mean to shift the curve accordingly.

#################################################################################################

peak_id = np.where(amp == np.max(amp))[0][0]

mean = pos_m[peak_id]

print mean

dif = pos_m - mean

t_d =(2.*dif)/(2.99*10**8)*10**(12.) # ps

print t_d

t = np.arange(0.,0.5,1/float(2*len(t_d)))        # picoseconds (ps)

################################

# Define constants

################################

c = 2.99*10**8                                  # m/s

alpha = np.pi                                   # rad

gamma = 200.0*10**(-15.)

lamb_o = 1160.00*10**(-9.)                      # m

omega_o = ((2.*np.pi*c)/lamb_o)*1.*10**(-12.)   # ps^-1

delta = np.pi/2.                                # rad

FWHM = 32.53*10**(-6)                           # m

t_p = ((FWHM)/c)*10**(12.)                      # ps

print t_p

E_norm = 1.                   # normalized

################################

# Define functions

################################

################################

# Electric field

################################

def E_o(x):

return E_norm*(cosh(1.76*x/t_p))**(-1.)

# Real part of electric field

def E_rin(x):

return (1./2.)*E_o(x)*cos(omega_o*x)

# Imaginary part of electric field

def E_iin(x):

return (1./2.)*E_o(x)*sin(omega_o*x)

# Total electric field

def E_(x):

return (1./2.)*E_o(x)*np.exp(-1j*omega_o*x)

################################

# Autocorrelation

################################

'''

def integrand(t,t_d):

return abs(E_rin(t))**2.*abs(E_rin(t - t_d))**2.

'''

def integrand(x,y):

return abs(E_(x))**2.*abs(E_(x - y))**2.

#integrand = abs(E_(t))**2.*abs(E_(t - t_d))**2.

def G(y):

return quad(integrand, -np.infty, np.infty, args=(y))

G_plot = []

for tau in t_d:

integral,error = G(tau)

G_plot.append(integral)

#fit data

params = curve_fit(fit,pos[174-100:174+100],amp[174-100:174+100],p0=[0.003,
8550,350]) #function, xdata, ydata, initial guess (from plot)

#parameters

[a,b,d] = params[0]

#error

[da,db,dd] = np.sqrt(np.diag(params[1]))

#################################

# Shaped electric field

#################################

alpha = np.pi

delta = np.pi/2.

gamma = 200.*10**(-15)       # sec

dt = (t[1]-t[0])*(1*10**(-12))

def phi(x):

return alpha*np.cos(gamma*(x - omega_o) - delta)

def M(x):

return np.exp(1j*phi(x))

def E_out(x):

E_in_w = fft(E_(x))

omega = fftfreq(len(x),dt)*2*np.pi

E_out_w = E_in_w*M(omega)

return ifft(E_out_w)

#################################

# Second autocorrelation

#################################

def integrand1(x,y):

return abs(E_out(x))**2.*abs(E_out(x - y))**2.

def G1(y):

return quad(integrand1, -np.infty, np.infty, args=(y))

G_plot_1 = []

for tau in t_d:

integral,error = G1(tau)

G_plot_1.append(integral)

#################################

# Plotting data

#################################

# Defining x and y minorlocator

xminorLocator = AutoMinorLocator()

yminorLocator = AutoMinorLocator()

# Setting minor ticks

ax = plt.subplot(111)

ax.xaxis.set_minor_locator(xminorLocator)

ax.yaxis.set_minor_locator(yminorLocator)

sample_pos = np.linspace(min(pos),max(pos),1000)

equation_label = str(round(a*1E3,1)) + r'$\times 10^{-3} \exp((x-$' +
str(round(b,1)) + r'$)/$' + str(round(d,1)) + r'$)$'

'''

subplot(2,1,1)

plot(pos,amp, 'r.', label='Data')

plot(sample_pos,fit(sample_pos,a,b,d),'k-', label='Guassian Fit')

title('Autocorrelation Intensity')

#xlabel('Stage Position (um)')

ylabel('Amplitude (V)')

text(8600,0.003,equation_label)

text(8800 , 0.0025 ,'Transform Limited Pulse',

horizontalalignment='center',

verticalalignment='center')

legend(loc='upper left')

subplot(2,1,2)

plot(pos, amp_auto, 'b-')

#title('Autocorrelation Intensity')

xlabel('Stage Position (um)')

ylabel('Amplitude (V)')

text(8800 , 0.0005 , r'$\Phi_M(\omega) = \alpha\cos(\gamma(\omega - \omega_0) - \delta)$',

horizontalalignment='center',

verticalalignment='center')

show()

'''

plot(t,E_(t))

title('Electric Field')

xlabel('Time (ps)')

ylabel('E (V/m)')

text(0.35, 0.35, r'$E_{out}(t)=E_{in}(t)=\frac{1}{2}E_0(t)e^{-i\omega_0t}$',

horizontalalignment='center',

verticalalignment='center')

'''

subplot(1,2,2)

plot(t,E_iin(t))

title('Electric Field')

xlabel('Time (sec)')

'''

show()

plot(t_d,G_plot)

title('Non-collinear Autocorrelation')

xlabel('Time (ps)')

ylabel('G(t_d)')

text(12.5, 0.8e-16, r'$G(t_d) = \int_{-\infty}^{\infty} |E(t)|^2 |E(t - t_d)|^2 dt$',

horizontalalignment='center',

verticalalignment='center')

show()

plot(t_d,G_plot_1)

xlabel('Time (ps)')

ylabel('G(t_d)')

text(12.5, 0.8e-16, r'$G(t_d) = \int_{-\infty}^{\infty} |E(t)|^2 |E(t - t_d)|^2 dt$',

horizontalalignment='center',

verticalalignment='center')

show()

plot(t,E_out(t))

title('Shaped Electric Field')

xlabel('Time (ps)')

ylabel('E_out')

show()

The stack trace is as follows:

IndexError                                Traceback (most recent call last)

<ipython-input-181-9d12acba3760> in <module>()

----> 1 G1(t_d[0])

<ipython-input-180-4eb6a4e577f9> in G1(y)

7 def G1(y):

8

----> 9     return quad(integrand1, -np.infty, np.infty, args=(y))

in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points,
weight, wvar, wopts, maxp1, limlst)

279         args = (args,)

280     if (weight is None):

--> 281         retval =

282     else:

283         retval =

in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)

345             return

346         else:

--> 347             return

348     else:

349         if infbounds != 0:

<ipython-input-180-4eb6a4e577f9> in integrand1(x, y)

1 def integrand1(x,y):

2

----> 3     return abs(E_out(x))**2.*abs(E_out(x - y))**2.

4

5

<ipython-input-144-8048dc460766> in E_out(x)

6

7 def E_out(x):

----> 8     E_in_w = fft(E_(x))

9     omega = fftfreq(len(x),dt)*2*np.pi

10     E_out_w = E_in_w*M(omega)

/Users/colinross/anaconda/lib/python2.7/site-packages/scipy/fftpack/basic.pyc
in fft(x, n, axis, overwrite_x)

253

254     if n is None:

--> 255         n = tmp.shape[axis]

256     elif n != tmp.shape[axis]:

257         tmp, copy_made = _fix_shape(tmp,n,axis)

IndexError: tuple index out of range

On Fri, Mar 13, 2015 at 3:27 PM, Danny Yoo <dyoo at hashcollision.org> wrote:

> On Fri, Mar 13, 2015 at 11:00 AM, Danny Yoo <dyoo at hashcollision.org>
> wrote:
> >> The error I am recieving is as follows:
> >>
> >> TypeError: only length-1 arrays can be converted to Python scalars
> >
> >
> > Hi Colin,
> >
> > Do you have a more informative "stack trace" of the entire error?
> > Providing this will help localize the problem.  As is, it's clear
> > there's a type error... somewhere... :P  Seeing the stack trace will
> > make things easier for us.
>
>
> Quick and dirty analysis.  If I had to guess, I'd look at the mix of
> math.erfc with numeric array values in the subexpression:
>
>     math.erfc((np.sqrt(2.)*x*1.E-3)/b)
>
>
> np.sqrt(...) returns a numpy array, if I'm not mistaken.  But
> math.erfc() probably can't deal with numpy values.  So there's
> definitely a problem there.
>
>
> See messages such as:
>
>
> which suggest that trying to use the standard library math functions
> on numpy arrays isn't going to work.
>
>
> Unfortunately, without explicit stack trace, I don't know if that's
> the *only* problem you're seeing.  That's why providing a good stack
> trace is so important in bug reports: there can be multiple causes for
> something to go wrong.  I'd rather make sure we've hit the one that's
> causing you grief.
>
>
>
> Anyway, remediation time.  Reading docs... ok, you could probably make
> a vectorized version of the erfc function:
>
>     vectorized_erfc = np.vectorize(math.erfc)
>
> and use this in favor of the non-vectorized version.  vectorize allows
> you to take regular functions and turn them into ones that work on
> numpy arrays:
>
>
> http://docs.scipy.org/doc/numpy/reference/generated/numpy.vectorize.html
>
>
> A better approach is probably to reuse the scipy.special.erfc function
> within scipy itself:
>
>
> http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.special.erfc.html
>
>
>
> If you have more questions, please feel free to ask.
>