Re: [Tutor] Reading in large .dat file
On Wed, Jun 29, 2016 at 4:21 PM, Peter Otten <__pete...@web.de> wrote: > Colin Ross wrote: > > > On Wed, Jun 29, 2016 at 2:41 PM, Peter Otten <__pete...@web.de> wrote: > > > >> Colin Ross wrote: > >> > >> > Good afternoon, > >> > > >> > I have a .dat file that is 214 rows by 65 columns that I would like to > >> > read into python and plot the 1st column versus all other columns. > >> > > >> > My code: > >> > > >> > > # > >> > > >> > import numpy as np > >> > import scipy > >> > import pylab as pl > >> > import matplotlib > >> > from matplotlib.ticker import ScalarFormatter, FormatStrFormatter > >> > import sys > >> > > >> > # Load in text from .dat file > >> > > >> > sed = np.loadtxt('spectra.dat', unpack = True) > >> > >> for flux in sed: > >> pl.plot(wavelength, flux) > >> > >> > pl.xscale('log') > >> > > >> > pl.show() > >> > > >> > > # > >> > > >> > This is fine if I want to write out a separate line for each of the 65 > >> > columns, but I would like to simplify the code by looping over the > >> > data. Can someone please help me formatting the loop correctly? > >> > > > > Thank you for the fast response! I have tried this and it outputs the > > attached image. > > This is a text-only mailing list, so the image didn't make it through. > > > It seems to only plot 1 curve, when I would have expected > > 64? > > You probably got a traceback (error information on the commandline). It is > always a good idea to include that. Example: > > $ python3 flux2.py > Traceback (most recent call last): > File "flux2.py", line 15, in > pl.plot(wavelength,flux) > NameError: name 'wavelength' is not defined > > The error indicates that there is no variable called wavelength. That's > because I didn't note that you use the first column to define that, sorry. > Here's a modified script that should work: > > import numpy as np > import pylab as pl > > sed = np.loadtxt('spectra.dat', unpack=True) > > wavelength = sed[0] > for flux in sed[1:]: > pl.plot(wavelength, flux) > > pl.xscale('log') > pl.show() > > The for loop now skips the first column. > Ah yes, this works perfect now. I had added the wavelength = sed[0], but forgot to skip it in the loop! Thank you! > > > ___ > Tutor maillist - Tutor@python.org > To unsubscribe or change subscription options: > https://mail.python.org/mailman/listinfo/tutor > ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
Re: [Tutor] Reading in large .dat file
On Wed, Jun 29, 2016 at 2:41 PM, Peter Otten <__pete...@web.de> wrote: > Colin Ross wrote: > > > Good afternoon, > > > > I have a .dat file that is 214 rows by 65 columns that I would like to > > read into python and plot the 1st column versus all other columns. > > > > My code: > > > > # > > > > import numpy as np > > import scipy > > import pylab as pl > > import matplotlib > > from matplotlib.ticker import ScalarFormatter, FormatStrFormatter > > import sys > > > > # Load in text from .dat file > > > > sed = np.loadtxt('spectra.dat', unpack = True) > > for flux in sed: > pl.plot(wavelength, flux) > > > pl.xscale('log') > > > > pl.show() > > > > # > > > > This is fine if I want to write out a separate line for each of the 65 > > columns, but I would like to simplify the code by looping over the data. > > Can someone please help me formatting the loop correctly? > Thank you for the fast response! I have tried this and it outputs the attached image. It seems to only plot 1 curve, when I would have expected 64? > > To understand why the above works (assuming it does) consider that for many > data types > > for item in data: > ... # use item > > is equivalent to > > for index in range(len(data)): > item = data[index] > ... # use item > > Unrolling the above for len(data) == 3 gives: > > item = data[0] > ... # use first item > item = data[1] > ... # use second item > item = data[2] > ... # use third item > > > ___ > Tutor maillist - Tutor@python.org > To unsubscribe or change subscription options: > https://mail.python.org/mailman/listinfo/tutor > ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] Reading in large .dat file
Good afternoon, I have a .dat file that is 214 rows by 65 columns that I would like to read into python and plot the 1st column versus all other columns. My code: # import numpy as np import scipy import pylab as pl import matplotlib from matplotlib.ticker import ScalarFormatter, FormatStrFormatter import sys # Load in text from .dat file sed = np.loadtxt('spectra.dat', unpack = True) # Access column data col1 = sed[0] col2 = sed[1] col3 = sed[2] col4 = sed[3] col5 = sed[4] wavelength = col1 flux_1 = col2 flux_2 = col3 flux_3 = col4 flux_4 = col5 pl.plot(wavelength,flux_1) pl.plot(wavelength,flux_2) pl.plot(wavelength,flux_3) pl.plot(wavelength,flux_4) pl.xscale('log') pl.show() # This is fine if I want to write out a separate line for each of the 65 columns, but I would like to simplify the code by looping over the data. Can someone please help me formatting the loop correctly? Thank you. Colin ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
Re: [Tutor] Extracting bits from an array
Thank you Michael, this is exactly what I was looking for! Clears things up on multiple fronts : ) On Fri, Apr 29, 2016 at 2:29 PM, Michael Selik <michael.se...@gmail.com> wrote: > > > > On Apr 29, 2016, at 1:15 PM, Colin Ross <colin.ross@gmail.com> > wrote: > > > > Hello, > > > > I have an array that takes on the following form: > > > > x = [1000,1001,1011,] > > > > The array elements are meant to be binary representation of integers. > > > > Goal: Access array elements and extract the first two bits. > > > > e.g. Final result would look something like this: > > > > x_new = [10,10,10,11] > > > > What I have tried: > > > > data_indices = range(4) # Set up array of values to loop over > > > > for idx in data_indices: > > f = x[idx] # Index into array of x values > > Instead of looping over a range of indices, you should loop over the data > itself. > > for number in x: > s = bin(number) > print s > > > f_idx = f[:2] # Extract first two elements > > You couldn't slice an integer. First convert to the binary representation > in string form. You can strip off the prefix if you just want the digits. > > s = bin(number).lstrip('0b') > > Then you can slice off the first two digits if you want. Remember, it's a > str of digits, not a number. > > > print f_idx > > > > I then receive the following error: > > > > IndexError: invalid index to scalar variable. > > > > Any help with accomplishing my outline dgoal would be greatly > appreciated. > > > > Thank you. > > > > Colin > > ___ > > Tutor maillist - Tutor@python.org > > To unsubscribe or change subscription options: > > https://mail.python.org/mailman/listinfo/tutor > ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] Extracting bits from an array
Hello, I have an array that takes on the following form: x = [1000,1001,1011,] The array elements are meant to be binary representation of integers. Goal: Access array elements and extract the first two bits. e.g. Final result would look something like this: x_new = [10,10,10,11] What I have tried: data_indices = range(4) # Set up array of values to loop over for idx in data_indices: f = x[idx] # Index into array of x values f_idx = f[:2] # Extract first two elements print f_idx I then receive the following error: IndexError: invalid index to scalar variable. Any help with accomplishing my outline dgoal would be greatly appreciated. Thank you. Colin ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] how to define constant width errorbars on log-log plot
Does anyone have an example code that shows how to plot errorbars with a constant line width for a dataset on a log log plot? Thank you. ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] Plotting asymmetric error bars for a single point in matplotlib
Hi all, Goal: To plot asymmetric x error bars for a single point using errorbar. I am interested in displaying the inter quartile range (IQR) for a data set. Code: import numpy as np import matplotlib.pyplot as plt y = 1.0 data = np.random.rand(100) median = np.median(data) upper_quartile = np.percentile(data, 75) lower_quartile = np.percentile(data, 25) IQR = upper_quartile - lower_quartile plt.errorbar(median, y, xerr=[lower_quartile ,upper_quartile], fmt='k--') plt.savefig('IQR.eps') plt.show() Error: Traceback (most recent call last): File IQR.py, line 15, in module plt.errorbar(median, y, xerr=[0.5,0.75], fmt='k--') File /usr/lib/pymodules/python2.7/matplotlib/pyplot.py, line 2251, in errorbar ret = ax.errorbar(x, y, yerr, xerr, fmt, ecolor, elinewidth, capsize, barsabove, lolims, uplims, xlolims, xuplims, **kwargs) File /usr/lib/pymodules/python2.7/matplotlib/axes.py, line 5327, in errorbar in cbook.safezip(x,xerr)] File /usr/lib/pymodules/python2.7/matplotlib/cbook.py, line 1294, in safezip raise ValueError(_safezip_msg % (Nx, i+1, len(arg))) ValueError: In safezip, len(args[0])=1 but len(args[1])=2 My understanding is that safezip zips together the x array with the specified upper and lower x limits? Thanks. Colin ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
Re: [Tutor] Shading Between Curves with Different Colour Over Specified X value Range
On Tue, Jul 28, 2015 at 11:03 AM, Oscar Benjamin oscar.j.benja...@gmail.com wrote: On Mon, 27 Jul 2015 at 20:53 Colin Ross colin.ross@gmail.com wrote: *Goal:* Shade between I_2 (curve 1) and I_3 (curve 2) with following conditions: - Green for 0 x 4 - Red for 4 x 12 *Code: * *Note: Code currently only attempting to shade green for 0 x 4 * import numpy as np import pylab from pylab import * import matplotlib.pyplot as plt import csv # Load data from .txt file with open('current_mirror_output_swing.csv', 'rb') as f: reader = csv.reader(f) your_list = list(reader) data = np.asarray(your_list) I_ref = np.asarray(data[1:,0]) I_1 = data[1:,1] I_2 = data[1:,2] I_3 = data[1:,3] # Create an array of x values to fill b/w curves with a certain color. Here you specify the X values for the fill shape as going from 0 to 4: X1 = np.linspace(0.,4.,len(I_3)) I_ref = I_ref.astype(float)*1000. I_1 = I_1.astype(float)*1000. I_2 = I_2.astype(float)*1000. I_3 = I_3.astype(float)*1000. # Plotting commands. Here you specify the X values for the line plots as being whatever I_ref is (I don't know what it is since I can't see your csv file): plot(I_ref, I_2, 'r-') plot(I_ref, I_3, 'b-') title('Current Mirror Output Swing') xlabel('$I_{ref}$ (mA)') ylabel('I (mA)') Try changing this line: plt.fill_between(X1, I_2, I_3, color = 'g', alpha = '0.5') to: plt.fill_between(I_ref, I_2, I_3, color = 'g', alpha = '0.5') and then the filled area should have the same X range as the lines. -- Oscar Thanks! This works now. ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] Arrow of constant size on log plot
I am attempting to draw an arrow that begins at a specified point on a logarithmic plot and then extends a certain distance in the -Y direction. This part is fine, however I would like to draw multiple arrows with the same size corresponding to different points on a plot. If I specify a certain value (as shown in the example below), the arrow size will change depending on where the point is on the plot, since the X and Y axis are logarithmic. - import numpy as np import pylab as pl x_1 = np.arange(10) y_1 = np.arange(10) x_1_avg = np.sum(x_1)/len(x_1) y_1_avg = np.sum(y_1)/len(y_1) x_2 = 3. y_2 = 3. pl.plot(x_1,y_1,'k') pl.arrow(x_1_avg,y_1_avg , 0.0, -0.5, fc='b', ec='b', head_width=0.1, head_length=0.1) pl.arrow(x_2,y_2, 0.0, -0.5, fc='r', ec='r', head_width=0.1, head_length=0.1 ) pl.yscale('log') pl.xscale('log') pl.show() - Any advice is appreciated. ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
Re: [Tutor] Shading Between Curves with Different Colour Over Specified X value Range
On Mon, Jul 27, 2015 at 5:01 PM, Mark Lawrence breamore...@yahoo.co.uk wrote: On 27/07/2015 19:47, Colin Ross wrote: *Goal:* Shade between I_2 (curve 1) and I_3 (curve 2) with following conditions: - Green for 0 x 4 - Red for 4 x 12 *Code: * *Note: Code currently only attempting to shade green for 0 x 4 * import numpy as np import pylab from pylab import * import matplotlib.pyplot as plt import csv # Load data from .txt file with open('current_mirror_output_swing.csv', 'rb') as f: reader = csv.reader(f) your_list = list(reader) data = np.asarray(your_list) I_ref = np.asarray(data[1:,0]) I_1 = data[1:,1] I_2 = data[1:,2] I_3 = data[1:,3] # Create an array of x values to fill b/w curves with a certain color. X1 = np.linspace(0.,4.,len(I_3)) I_ref = I_ref.astype(float)*1000. I_1 = I_1.astype(float)*1000. I_2 = I_2.astype(float)*1000. I_3 = I_3.astype(float)*1000. # Plotting commands. plot(I_ref, I_2, 'r-') plot(I_ref, I_3, 'b-') title('Current Mirror Output Swing') xlabel('$I_{ref}$ (mA)') ylabel('I (mA)') plt.fill_between(X1, I_2, I_3, color = 'g', alpha = '0.5') plt.legend(['$I_{ref}$', '$I_{out}$'], loc='upper left') plt.grid() show() *Issue: * See attached figure. Thank you. There is no attachment to see, sorry :( My apologies. SHould be there now! One thing to note about the following lines. from pylab import * import matplotlib.pyplot as plt The first was designed to make matplotlib easy to use interactively, especially in iPython, the second in a script. IIRC the former is deprecated so I suggest you stick with the latter. Great, thank you! -- My fellow Pythonistas, ask not what our language can do for you, ask what you can do for our language. Mark Lawrence ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] Shading Between Curves with Different Colour Over Specified X value Range
*Goal:* Shade between I_2 (curve 1) and I_3 (curve 2) with following conditions: - Green for 0 x 4 - Red for 4 x 12 *Code: * *Note: Code currently only attempting to shade green for 0 x 4 * import numpy as np import pylab from pylab import * import matplotlib.pyplot as plt import csv # Load data from .txt file with open('current_mirror_output_swing.csv', 'rb') as f: reader = csv.reader(f) your_list = list(reader) data = np.asarray(your_list) I_ref = np.asarray(data[1:,0]) I_1 = data[1:,1] I_2 = data[1:,2] I_3 = data[1:,3] # Create an array of x values to fill b/w curves with a certain color. X1 = np.linspace(0.,4.,len(I_3)) I_ref = I_ref.astype(float)*1000. I_1 = I_1.astype(float)*1000. I_2 = I_2.astype(float)*1000. I_3 = I_3.astype(float)*1000. # Plotting commands. plot(I_ref, I_2, 'r-') plot(I_ref, I_3, 'b-') title('Current Mirror Output Swing') xlabel('$I_{ref}$ (mA)') ylabel('I (mA)') plt.fill_between(X1, I_2, I_3, color = 'g', alpha = '0.5') plt.legend(['$I_{ref}$', '$I_{out}$'], loc='upper left') plt.grid() show() *Issue: * See attached figure. Thank you. ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] Python serial interface
Hi all, This is a very general question, but I was wondering if anyone has experience using python to interface with a serial port? If so, can you please forward any useful resources? Thanks! Colin ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
Re: [Tutor] Python serial interface
I am using the following controller: http://www.aerotech.com/product-catalog/drives-and-drive-racks/ensemble-mp.aspx Which does not specifically list python as one of the accepted languages, but I guess this does not mean it is not possible. Colin On Wed, Apr 1, 2015 at 12:58 PM, Colin Ross colin.ross@gmail.com wrote: Thank you Francois, this gives me a lot to think about! I really appreciate your feedback. Colin On Wed, Apr 1, 2015 at 12:50 PM, Francois Dion francois.d...@gmail.com wrote: On Wed, Apr 1, 2015 at 11:01 AM, Colin Ross colin.ross@gmail.com wrote: Hi Francois, Thank you for the fast reply! I am looking to control a brushless servo motor ( http://www.aerotech.com/product-catalog/motors/rotary-motors/bms-series.aspx) that drives a rotary stage. These motors are not controlled by serial, you'll need a brushless controller. In turn the controller might be or not serial, or ip or something else altogether (PWM is typical). They do have an rs-422 connection, but that is for a feedback signal. So you are building a feedback control system. If you google that and python you'll find several useful references, including a design and analysis module (but not control): https://pypi.python.org/pypi/control/0.6.6 http://sourceforge.net/projects/python-control/ Nothing out of the box for your application, obviously. Also, the motor control itself will depend on the microprocessor hardware you are using and the motor controller. Assuming PWM, the raspberry pi has a software PWM module in Python. Same with micropython. Anything else is pretty much DIY. As for the overall concept of control systems, Chapter 9 of Real World Instrumentation ( http://shop.oreilly.com/product/9780596809577.do ) will give you an overview of what is involved (with some python example, but not directly applicable to your system). Francois -- raspberry-python.blogspot.com - www.pyptug.org - www.3DFutureTech.info http://www.3dfuturetech.info/ - @f_dion ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
Re: [Tutor] Python serial interface
Hi Francois, Thank you for the fast reply! I am looking to control a brushless servo motor ( http://www.aerotech.com/product-catalog/motors/rotary-motors/bms-series.aspx) that drives a rotary stage. Colin On Wed, Apr 1, 2015 at 11:53 AM, Francois Dion francois.d...@gmail.com wrote: Pyserial is python 2.x and 3.x compatible. It is very widely used and is stable. http://pyserial.sourceforge.net/ What is your application? Sometimes you can use a higher level module that makes use of pyserial. Francois -- raspberry-python.blogspot.com - www.pyptug.org - www.3DFutureTech.info - @f_dion ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
Re: [Tutor] Python serial interface
Thank you Francois, this gives me a lot to think about! I really appreciate your feedback. Colin On Wed, Apr 1, 2015 at 12:50 PM, Francois Dion francois.d...@gmail.com wrote: On Wed, Apr 1, 2015 at 11:01 AM, Colin Ross colin.ross@gmail.com wrote: Hi Francois, Thank you for the fast reply! I am looking to control a brushless servo motor ( http://www.aerotech.com/product-catalog/motors/rotary-motors/bms-series.aspx) that drives a rotary stage. These motors are not controlled by serial, you'll need a brushless controller. In turn the controller might be or not serial, or ip or something else altogether (PWM is typical). They do have an rs-422 connection, but that is for a feedback signal. So you are building a feedback control system. If you google that and python you'll find several useful references, including a design and analysis module (but not control): https://pypi.python.org/pypi/control/0.6.6 http://sourceforge.net/projects/python-control/ Nothing out of the box for your application, obviously. Also, the motor control itself will depend on the microprocessor hardware you are using and the motor controller. Assuming PWM, the raspberry pi has a software PWM module in Python. Same with micropython. Anything else is pretty much DIY. As for the overall concept of control systems, Chapter 9 of Real World Instrumentation ( http://shop.oreilly.com/product/9780596809577.do ) will give you an overview of what is involved (with some python example, but not directly applicable to your system). Francois -- raspberry-python.blogspot.com - www.pyptug.org - www.3DFutureTech.info http://www.3dfuturetech.info/ - @f_dion ___ 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
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
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
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
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)
[Tutor] Fitting data to error function
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
[Tutor] selecting data from a list
Hi all, I am attempting to isolate a certain subset of data from the list a and then turn it into any array. To do so, I have written the following code: import numpy as np a = [0,1,2,3,4,5,6,7,8,9,10] b = [10,20,30,40,50,60,70,80,90,100,110] for a in range(len(a)): if a 5: print a a_1 = np.array(a) print a_1 The output is as follows: 6 7 8 9 10 10 As you can see, when I attempt to turn the list of numbers 6 through 10 into an array I am left with it only printing out 10... The desired result is: [6,7,8,9,10} Any guidance would be greatly appreciated. Thank you. Colin ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] Plotting subsets of data
Good afternoon, I am using the following to code to plot the output from an optical encoder: 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 sys # Defining x and y minorlocator xminorLocator=AutoMinorLocator() yminorLocator=AutoMinorLocator() # Load data from .txt file data = np.loadtxt('2014_12_04-16_30_03.txt',skiprows = 0 ,usecols = (0,1)) print \n Chopper Test for X-SPEC Prototype print \n Time (sec) Pos (deg) print data # Place column data in array to be plotted time = np.array(data[:,0]) print Time (sec): #print time pos = np.array(data[:,1]) print Position (deg): print pos # Setting minor ticks ax = plt.subplot(111) ax.xaxis.set_minor_locator(xminorLocator) ax.yaxis.set_minor_locator(yminorLocator) #Plotting commands. subplot(2,1,1) plot(time, pos, 'ro') title('Encoder Output') ylabel('Pos (deg)') subplot(2,1,2) plot(t_small, pos, 'ro') xlabel('Time (sec)') ylabel('Pos (deg)') show() The desired result is a square wave, but this is not readily available from the plot (see attached). For the subplot(2,1,2) I would like to plot the output over a 5 second interval so that the behaviour becomes more evident. Can someone please advise me as to the easiest way isolate the first 5 seconds of data? Thank you. ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
Re: [Tutor] Plotting subsets of data
Perfect, thank you Danny! On Mon, Dec 8, 2014 at 3:49 PM, Danny Yoo d...@hashcollision.org wrote: On Mon, Dec 8, 2014 at 7:48 AM, Colin Ross colin.ross@gmail.com wrote: Good afternoon, I am using the following to code to plot the output from an optical encoder: Hi Colin, Matplotlib is a third-party library, so you may also consider asking the matplotlib folks. From a brief look at: http://matplotlib.org/1.4.2/users/pyplot_tutorial.html#working-with-multiple-figures-and-axes it appears that you can override the default axis, and specify xmin, ymin, xmax, and ymax values. http://matplotlib.org/1.4.2/api/pyplot_api.html#matplotlib.pyplot.axis For your particular case, you may want to just limit your x axis, in which case xlim() might be appropriate. http://matplotlib.org/1.4.2/api/pyplot_api.html#matplotlib.pyplot.xlim If all else fails, just filter your data before submitting it to the grapher. The program loads data here, using loadtxt (http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html): data = np.loadtxt('2014_12_04-16_30_03.txt',skiprows = 0 ,usecols = (0,1)) and it's just a numpy array: you can manipulate numpy arrays. See: http://stackoverflow.com/questions/26154711/filter-rows-of-a-numpy-array as an example of an approach. ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
[Tutor] Error message received when running “from scipy import interpolate”
I am attempting to run the following in python: from scipy import interpolate I receive this error message: Traceback (most recent call last): File stdin, line 1, in module File /Library/Python/2.7/site-packages/scipy-0.11.0.dev_0496569_20111005-py2.7-macosx-10.7-x86_64.egg/scipy/interpolate/__init__.py, line 156, in module from ndgriddata import * File /Library/Python/2.7/site-packages/scipy-0.11.0.dev_0496569_20111005-py2.7-macosx-10.7-x86_64.egg/scipy/interpolate/ndgriddata.py, line 9, in module from interpnd import LinearNDInterpolator, NDInterpolatorBase, \ File numpy.pxd, line 172, in init interpnd (scipy/interpolate/interpnd.c:7696)ValueError: numpy.ndarray has the wrong size, try recompiling I am currently running Mac OS X Lion 10.7.5 and Python 2.7.1. with MacPorts installed. I am attempting to run the code on a different computer than it was created on. After doing some reading my understanding is that this error may be a result of the ABI not being forward compatible. To try and solve the issue I updated my MacPorts and then ran ($ sudo port upgrade outdated) to upgrade the installed ports. However, the same error continues to appear when I try and run the code. ___ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor