[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
from scipy.integrate import quad,fixed_quad
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)
title('Masked')
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))
/Users/colinross/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc
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 =
_quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
282 else:
283 retval =
_quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)
/Users/colinross/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc
in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
345 return
_quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
346 else:
--> 347 return
_quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
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:
>
> https://groups.google.com/forum/#!topic/comp.lang.python/9E4HX4AES-M
>
> 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.
>
More information about the Tutor
mailing list