Re: [Tutor] Fitting data to error function

2015-03-16 Thread Oscar Benjamin 
On 16 March 2015 at 22:22, Colin Ross colin.ross@gmail.com wrote:

 Yes, thank you, they were suppose to both be E_out.

Hi Colin,

I'm not sure if that means that your problem is fixed or not but I
thought I would point something out that helps in fixing this kind of
problem.

You're using ipython which has an excellent debugger. So let's say I
have the following buggy (and not very useful) program called tmp.py:

from numpy.fft import fft

def g(z):
return str(z)

def f(x):
return fft(g(x))

f(1)

I can run the program under ipython by typing run tmp.py:

$ ipython
Python 2.7.6 (default, Mar 22 2014, 22:59:56)
Type copyright, credits or license for more information.

IPython 1.2.1 -- An enhanced Interactive Python.
? - Introduction and overview of IPython's features.
%quickref - Quick reference.
help  - Python's own help system.
object?   - Details about 'object', use 'object??' for extra details.

In [1]: run tmp.py
---
IndexErrorTraceback (most recent call last)
/usr/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in
execfile(fname, *where)
202 else:
203 filename = fname
-- 204 __builtin__.execfile(filename, *where)

/stuff/oscar/work/current/tmp/ipdb/tmp.py in module()
  7 return fft(g(x))
  8
 9 f(1)

/stuff/oscar/work/current/tmp/ipdb/tmp.py in f(x)
  5
  6 def f(x):
 7 return fft(g(x))
  8
  9 f(1)

/usr/lib/python2.7/dist-packages/numpy/fft/fftpack.pyc in fft(a, n, axis)
172 
173
-- 174 return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf,
_fft_cache)
175
176

/usr/lib/python2.7/dist-packages/numpy/fft/fftpack.pyc in _raw_fft(a,
n, axis, init_function, work_function, fft_cache)
 48
 49 if n is None:
--- 50 n = a.shape[axis]
 51
 52 if n  1:

IndexError: tuple index out of range


That's great because the information that I see shows me where the
error occurs - if I know how to read it. In this case I look up until
I find a line that is part of my own code i.e. line 7 of tmp.py.
That's where I passed something into the numpy fft function which lead
to an error.

At this point I can try to understand what is the object that is being
passed into fft by looking at my code. Or I can use the ipython
debugger to manually inspect the object. Turn the debugger on by
typing pdb. Then run your script again:

In [2]: pdb
Automatic pdb calling has been turned ON

In [3]: run tmp.py
---
IndexErrorTraceback (most recent call last)
/usr/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in
execfile(fname, *where)
202 else:
203 filename = fname
-- 204 __builtin__.execfile(filename, *where)

/stuff/oscar/work/current/tmp/ipdb/tmp.py in module()
  7 return fft(g(x))
  8
 9 f(1)

/stuff/oscar/work/current/tmp/ipdb/tmp.py in f(x)
  5
  6 def f(x):
 7 return fft(g(x))
  8
  9 f(1)

/usr/lib/python2.7/dist-packages/numpy/fft/fftpack.pyc in fft(a, n, axis)
172 
173
-- 174 return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf,
_fft_cache)
175
176

/usr/lib/python2.7/dist-packages/numpy/fft/fftpack.pyc in _raw_fft(a,
n, axis, init_function, work_function, fft_cache)
 48
 49 if n is None:
--- 50 n = a.shape[axis]
 51
 52 if n  1:

IndexError: tuple index out of range
 /usr/lib/python2.7/dist-packages/numpy/fft/fftpack.py(50)_raw_fft()
 49 if n is None:
--- 50 n = a.shape[axis]
 51

ipdb u
 /usr/lib/python2.7/dist-packages/numpy/fft/fftpack.py(174)fft()
173
-- 174 return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf,
_fft_cache)
175

ipdb u
 /stuff/oscar/work/current/tmp/ipdb/tmp.py(7)f()
  6 def f(x):
 7 return fft(g(x))
  8

ipdb p x
1
ipdb p g(x)
'1'
ipdb p type(g(x))
type 'str'

The debugger with the ipdb prompt freezes the script at the point
where the error occurs and enables me to manually inspect what is
going on. It started in the numpy code so I typed u twice to go up
two frames to where my function called the function that called the
function that generated the exception. Here I can query the object x
(typing p for print) and test what g(x) returns. Now the problem
is clear to me: I'm passing a string into fft when it is expecting an
array of numbers.

The python interpreter comes with a slightly less convenient debugger
called pdb that you can use in the same way if you run your script
from the terminal like so:

$ python -m pdb tmp.py
 /stuff/oscar/work/current/tmp/ipdb/tmp.py(1)module()
- from numpy.fft import fft
(Pdb) c
Traceback (most recent call last):
  File /usr/lib/python2.7/pdb.py, line 1314, in

Re: [Tutor] Fitting data to error function

2015-03-16 Thread Colin Ross 
What I am trying to do is calculate the non-colinear autocorrelation:

G(t_d) = \int_{-\infty}^{+\infty} |E(t)|^2 * |E(t - t_d)|^2 dt

So I need to loop through an array of t_d values (len = 376) and calculate
G(t_d) for as many t values as possible to eliminate sampling issues.

Colin

On Mon, Mar 16, 2015 at 6:50 PM, Colin Ross colin.ross@gmail.com
wrote:

 HI Danny,

 Here is a simplified version:

 

 import numpy as np

 import pylab

 from pylab import *

 import matplotlib.pyplot as plt

 import scipy

 from scipy.integrate import quad

 from scipy.fftpack import fft, ifft, fftfreq

 ##

 # 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))

 

 # Create data arrays

 

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

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

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

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

 #

 # 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]

 dif = pos_m - mean

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


 t = np.arange(-0.5,0.5,0.0001)# 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


 E_norm = 1. #
 normalized

 #

 # Define functions

 #


 

 # Input electric field

 

 def E_o(x):

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

 def E_(x):

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

 #

 # Shaped electric field

 #

 alpha = np.pi

 delta = np.pi/2.

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

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

 def phi(x):

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

 def M(x):

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

 ##

 # Write E(w) as E(t) using FFT

 ##

 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)

 def integrand(x,y):

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

 def G(y):

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

 auto_corr = []

 for tau in t_d:

 integral,error = G(tau)

 auto_corr.append(integral)

 # Plotting

 plt.plot(t_d,auto_corr)

 plt.show()


 -

 ---

 IndexErrorTraceback (most recent call last)

 /Users/colinross/Documents/Optics/Programming/five.py in module()

 107

 108 for tau in t_d:

 -- 109 integral,error = G(tau)

 110 auto_corr.append(integral)

 111


 /Users/colinross/Documents/Optics/Programming/five.py in G(y)

 102

 103 def G(y):

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

 105

 106 auto_corr = []


 /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)

Re: [Tutor] Fitting data to error function

2015-03-16 Thread Colin Ross 
Yes, thank you, they were suppose to both be E_out.

And to answer your last question, I do not. Can you please explain?



On Mon, Mar 16, 2015 at 7:19 PM, Danny Yoo d...@hashcollision.org wrote:

 On Mon, Mar 16, 2015 at 2:55 PM, Colin Ross colin.ross@gmail.com
 wrote:
  What I am trying to do is calculate the non-colinear autocorrelation:
 
  G(t_d) = \int_{-\infty}^{+\infty} |E(t)|^2 * |E(t - t_d)|^2 dt
 
  So I need to loop through an array of t_d values (len = 376) and
 calculate
  G(t_d) for as many t values as possible to eliminate sampling issues.


 Ok.  But you're using the term E(t) and E(t-t_d) in your
 LaTeX-ified equation in such a way that it sounds like 'E' is context
 sensitive.

 Look at the Python definition of integrand() again:

 100 def integrand(x,y):
 -- 101 return abs(E_out(x))**2.*abs(E_(x - y))**2.


 and note that there are *two* distinct functions here being used:

 E_out
 E_

 In contrast, in your mathematics, it looks like these should be the
 *same* E function, so the fact that this is *different* is worth
 consideration.  I have to assume that the mathematics has some context
 sensitivity based on the argument type that you haven't quite
 explained here.


 Also, I need to ask: do you know what is meant by the term unit
 test?  Because this doesn't seem to have been addressed yet, so I
 need to double check.

___
Tutor maillist  -  Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Fitting data to error function

2015-03-16 Thread Colin Ross 
HI Danny,

Here is a simplified version:



import numpy as np

import pylab

from pylab import *

import matplotlib.pyplot as plt

import scipy

from scipy.integrate import quad

from scipy.fftpack import fft, ifft, fftfreq

##

# 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))



# Create data arrays



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

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

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

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

#

# 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]

dif = pos_m - mean

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


t = np.arange(-0.5,0.5,0.0001)# 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


E_norm = 1. #
normalized

#

# Define functions

#




# Input electric field



def E_o(x):

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

def E_(x):

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

#

# Shaped electric field

#

alpha = np.pi

delta = np.pi/2.

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

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

def phi(x):

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

def M(x):

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

##

# Write E(w) as E(t) using FFT

##

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)

def integrand(x,y):

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

def G(y):

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

auto_corr = []

for tau in t_d:

integral,error = G(tau)

auto_corr.append(integral)

# Plotting

plt.plot(t_d,auto_corr)

plt.show()


-

---

IndexErrorTraceback (most recent call last)

/Users/colinross/Documents/Optics/Programming/five.py in module()

107

108 for tau in t_d:

-- 109 integral,error = G(tau)

110 auto_corr.append(integral)

111


/Users/colinross/Documents/Optics/Programming/five.py in G(y)

102

103 def G(y):

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

105

106 auto_corr = []


/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:


/Users/colinross/Documents/Optics/Programming/five.py in integrand(x, y)

 99

100 def integrand(x,y):

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

102

103 def G(y):


/Users/colinross/Documents/Optics/Programming/five.py in E_out(x)

 93

 94 def E_out(x):

--- 95

Re: [Tutor] Fitting data to error function

2015-03-16 Thread Danny Yoo 
On Mon, Mar 16, 2015 at 2:55 PM, Colin Ross colin.ross@gmail.com wrote:
 What I am trying to do is calculate the non-colinear autocorrelation:

 G(t_d) = \int_{-\infty}^{+\infty} |E(t)|^2 * |E(t - t_d)|^2 dt

 So I need to loop through an array of t_d values (len = 376) and calculate
 G(t_d) for as many t values as possible to eliminate sampling issues.


Ok.  But you're using the term E(t) and E(t-t_d) in your
LaTeX-ified equation in such a way that it sounds like 'E' is context
sensitive.

Look at the Python definition of integrand() again:

100 def integrand(x,y):
-- 101 return abs(E_out(x))**2.*abs(E_(x - y))**2.


and note that there are *two* distinct functions here being used:

E_out
E_

In contrast, in your mathematics, it looks like these should be the
*same* E function, so the fact that this is *different* is worth
consideration.  I have to assume that the mathematics has some context
sensitivity based on the argument type that you haven't quite
explained here.


Also, I need to ask: do you know what is meant by the term unit
test?  Because this doesn't seem to have been addressed yet, so I
need to double check.
___
Tutor maillist  -  Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Fitting data to error function

2015-03-16 Thread Danny Yoo 
On Mon, Mar 16, 2015 at 3:22 PM, Colin Ross colin.ross@gmail.com wrote:
 Yes, thank you, they were suppose to both be E_out.

 And to answer your last question, I do not. Can you please explain?

The article:

http://en.wikipedia.org/wiki/Unit_testing

may help.


As a brief intro: we treat each interesting function as a
independently-testable unit, and write code that exercises each
function against a series of input-expected output pairs.

E.g. toy example: if we were writing a function to compute hypotenuse:

##
def h(x, y):
return sqrt(sq(x) + sq(y))

def sqrt(x):
return x**0.5

def sq(x):
return x ** 2
##

then we can unit test the individual functions.



Here is a complete toy example of the idea:

###
import unittest

def hypo(x, y):
return sqrt(sq(x) + sq(y))

def sqrt(x):
return x**0.5

def sq(x):
return x ** 2

class HypoTestCases(unittest.TestCase):
   def testSqrt(self):
self.assertEqual(2, sqrt(4))
self.assertEqual(0, sqrt(0))

   def testSq(self):
self.assertEqual(1, sq(1))
self.assertEqual(4, sq(2))
self.assertEqual(9, sq(3))

   def testHypo(self):
self.assertEqual(5, hypo(3, 4))

if __name__ == '__main__':
unittest.main()




When we run this program, the unit tests will execute and tell us if
they see any deviation between expectations and reality.


$ python hypo.py
...
--
Ran 3 tests in 0.000s

OK




Now let us intentionally damage one of the functions.  Change the
definition of sq(x): to:

#
def sq(x):
return x * 2
#

and execute the unit tests again.  We'll now see something very interesting:


##
$ python hypo.py
FF.
==
FAIL: testHypo (__main__.HypoTestCases)
--
Traceback (most recent call last):
  File hypo.py, line 23, in testHypo
self.assertEqual(5, hypo(3, 4))
AssertionError: 5 != 3.7416573867739413

==
FAIL: testSq (__main__.HypoTestCases)
--
Traceback (most recent call last):
  File hypo.py, line 18, in testSq
self.assertEqual(1, sq(1))
AssertionError: 1 != 2

--
Ran 3 tests in 0.000s

FAILED (failures=2)
##


Our unit test code will automatically report if the actual output of
a function mismatches against the expected output.  Flaws in auxiliary
functions can be reported, which means that if we do see problems, we
can more easily scope what the issue is.  Note: from the unit test
error output, we know that sqrt() is probably ok, but hypo() and sq()
are probably broken.  And since hypo()'s definition depends on sq(),
it naturally leads us to investigate sq() first.


It is the use of tests on individual units of our code that lets us
determine a fault without having to look at the entire aggregate
program.  Contrast this with the situation you're in at the moment:
you're currently stuck trying to debug the *entire program* when
something goes wrong.  Unfortunately, this strategy doesn't scale when
programs get too large.  An effective mechanism that programmers use
is to unit test the individual functions.
___
Tutor maillist  -  Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Fitting data to error function

2015-03-15 Thread Colin Ross 
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)

Re: [Tutor] Fitting data to error function

2015-03-15 Thread Danny Yoo 
 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.


At this point, the program is large enough that we need to apply
alternative techniques besides reading the entire program and trying
to keep it in our head at once.  One approach that is often successful
is to check the _user functions_ in isolation, because they should be
separable and independently testable.

You've written a few functions in your program.  In particular:

##
def fit(x,a,b,c): ...
def E_o(x): ...
def E_rin(x): ...
def E_iin(x): ...
def E_(x): ...
def integrand(t,t_d): ...
def integrand(x,y): ...## why is this being defined twice?!
def G(y): ...
def phi(x): ...
def M(x): ...
def E_out(x): ...
def integrand1(x,y): ...
def G1(y): ...
##

A function-based perspective will ask the following questions:

1.  What is the documented purpose of each of these functions?

2.  What are the expected inputs and output types for each of
these functions?

3.  Which of these functions have unit tests to spot-check their quality?

I can say for certain, without looking at any significant code, that
something is strange with regards to integrand.  Your program has
defined this twice, and that's almost certainly not right.

The strategy here is to pinpoint where the problem is: we can't
possibly say that the whole program is broken, so we try to scope the
problem down to a manageable size.

You mention that the first auto-correlation is producing good results.
What functions are involved in the computation of the first
auto-correlation?


Another technique to try to scope is to look at the stack trace and
the involved functions.  The stack trace says that the following
functions are involved at the point of the program breakage:


  G1
  quad
  integrand1
  E_out
  E_

In particular, note the tail end of the stack trace:



   7 def E_out(x):
  8 E_in_w = fft(E_(x))


It is on the call to:

fft(E_(x))

that things break with the IndexError.  Note that fft() is provided by
the scipy library.  For the purposes of debugging, let us assume that
all external third-party libraries are bug-free.

What does fft expect to receive as an argument?  We can read the following:


http://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.fft.html#scipy.fftpack.fft

Since fft is erroring out: there's only one possibility: E_(x) is not
providing a value that's appropriate to fft().

Why would it do that?  Two possibilities:

1.  E_ is buggy and needs investigation.

2.  E_ is fine, but the inputted value of x is one that E_x is not
defined to handle.


And so we start the investigation by considering those two
possibilities.  Start with #1.

Do you know if E_ is fine?  What is it supposed to do?  What is the
shape of its input?  What is the shape of its output?  Is its output
something that fft() is designed to handle?
___
Tutor maillist  -  Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Fitting data to error function

2015-03-15 Thread Danny Yoo 
 What does fft expect to receive as an argument?  We can read the following:

 
 http://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.fft.html#scipy.fftpack.fft

 Since fft is erroring out: there's only one possibility: E_(x) is not
 providing a value that's appropriate to fft().

 Why would it do that?  Two possibilities:

 1.  E_ is buggy and needs investigation.

 2.  E_ is fine, but the inputted value of x is one that E_x is not
 defined to handle.


 And so we start the investigation by considering those two
 possibilities.  Start with #1.

 Do you know if E_ is fine?  What is it supposed to do?  What is the
 shape of its input?  What is the shape of its output?  Is its output
 something that fft() is designed to handle?

Just to add: this is not to say that E_() is buggy.  We just have to
start investigation somewhere, and it seemed as good a place as any,
since I don't know if _any_ of the functions are behaving.  :P

This is very much the reason why programmers like to do unit testing:
we want to know what functions at least do something that we know is
useful.  We know all too well that whole programs break, and we want
to have some confidence on what components of our program are likely
ok.

If #2 is the actual issue, then the question becomes: why is the
program producing an invalid input 'x'?  And that's where we need to
start reasoning backwards, to discover how that value was constructed.
___
Tutor maillist  -  Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Fitting data to error function

2015-03-13 Thread Danny Yoo 
On Fri, Mar 13, 2015 at 11:00 AM, Danny Yoo d...@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.
___
Tutor maillist  -  Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor

[Tutor] Fitting data to error function

2015-03-13 Thread Colin Ross 
Hi all,

I am attempting to optimize the parameters I_0 and w_0 in the function
(func(x,I_0,w_0) shown below) to fit a set of data. However when I run this
code I recieve the error shown below. Any suggestions would be greatly
appreciated.

Code:

import numpy as np

import math

from math import *

# data

xdata = np.array([0.1,1.0,2.0,3.0,4.0,5.0])

ydata = np.array([0.1,0.9,2.2,2.8,3.9,5.1])

# Initial guess.

x0 = np.array([0.0, 0.0])


# Data errors:

sigma = np.array([1.0,1.0,1.0,1.0,1.0,1.0])

# Error function

def func(x,a,b):

return ((a)/2.)*(math.erfc((np.sqrt(2.)*x*1.E-3)/b))


# Optimization

import scipy.optimize as optimization

print optimization.curve_fit(func, xdata, ydata, x0, sigma)

The error I am recieving is as follows:

TypeError: only length-1 arrays can be converted to Python scalars
___
Tutor maillist  -  Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Fitting data to error function

2015-03-13 Thread Danny Yoo 
 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.
___
Tutor maillist  -  Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor