Interactive data analysis with iPython

Many users are familiar with matlab, which is not supported by the LCLS offline group. This is a description of how you can use interactive python to analyze your data in a very similar way.

Many thanks to the XPP crew, Henrik Lemke in particular, for providing guidance, code examples (in matlab) and feedback in writing the python scripts presented here.

The components

• ipython
  ... an interactive python shell
  External links:
  • http://ipython.github.com/ipython-doc/dev/interactive/index.html
  • http://ipython.scipy.org/moin/Cookbook
  • http://pages.physics.cornell.edu/~myers/teaching/ComputationalMethods/python/ipython.html

• matplotlib pyplot / pylab
  ... an open-source plotting tool

• numpy
  ... numeric python
  External links:
  • http://www.scipy.org/NumPy_for_Matlab_Users

The Matlab-like tools / library provided by PCDS

• pymatlab.py
  ... a module implementing in python some of the tools written by Henrik/XPP for matlab
  
  Starting iPython

  [ofte@psana0XXX myrelease]$ ipython
  Python 2.4.3 (#1, Nov  3 2010, 12:52:40)
  Type "copyright", "credits" or "license" for more information.
  
  IPython 0.9.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'. ?object also works, ?? prints more.
  

  
  Loading the library. Normally 'import pymatlab' would be recommended, but if you do 'from pymatlab import *', all the functions defined in this module gets loaded in the current namespace, and you can see them in your workspace. This might be easier for interactive work.

  In [1]: from pymatlab import *
  Pretend this is matlab
  
  

  
  who gives you a short list of workspace contents

  In [2]: who
  H5getobjnames   ScanInput       ScanOutput      filtvec findmovingmotor 
  getSTDMEANfrac_from_startpoint  get_filter get_limits       get_limits_automatic
  get_limits_channelhist  get_limits_correlation  get_limits_corrfrac     
  h5py    np      plt     rdXPPdata       runexpNO2fina       scan       scaninput
  

  
  whos gives you a longer list of workspace contents

  In [3]: whos
  Variable                         Type          Data/Info
  --------------------------------------------------------
  H5getobjnames                    function      <function H5getobjnames at 0x2b57de8>
  ScanInput                        type          <class 'pymatlab.ScanInput'>
  ScanOutput                       type          <class 'pymatlab.ScanOutput'>
  filtvec                          function      <function filtvec at 0x2b57f50>
  findmovingmotor                  function      <function findmovingmotor at 0x2b57d70>
  getSTDMEANfrac_from_startpoint   function      <function getSTDMEANfrac_<...>_startpoint at 0x2b581b8>
  get_filter                       function      <function get_filter at 0x2b57ed8>
  get_limits                       function      <function get_limits at 0x2b58050>
  get_limits_automatic             function      <function get_limits_automatic at 0x2b58230>
  get_limits_channelhist           function      <function get_limits_channelhist at 0x2b582a8>
  get_limits_correlation           function      <function get_limits_correlation at 0x2b580c8>
  get_limits_corrfrac              function      <function get_limits_corrfrac at 0x2b58140>
  h5py                             module        <module 'h5py' from '/reg<...>ython/h5py/__init__.pyc'>
  np                               module        <module 'numpy' from '/re<...>thon/numpy/__init__.pyc'>
  plt                              module        <module 'matplotlib.pyplo<...>n/matplotlib/pyplot.pyc'>
  rdXPPdata                        function      <function rdXPPdata at 0x2b57c80>
  runexpNO2fina                    function      <function runexpNO2fina at 0x2b57e60>
  scan                             ScanOutput    <pymatlab.ScanOutput object at 0x2b60536bee90>
  scaninput                        ScanInput     <pymatlab.ScanInput object at 0x2b60536b4e90>
  
  

Examples

1) Select limits from graphical input and plot filtered IPIMB

Here's a log from a session that produces a loglog plot (blue dots) of two IPIMB channels, selects limits from graphical inpu (mouse click),
draws the selected events with red dots.

In [1]: from pymatlab import *
Pretend this is matlab

In [2]: whos
Variable                         Type        Data/Info
------------------------------------------------------
H5getobjnames                    function    <function H5getobjnames at 0x1379fde8>
ScanInput                        type        <class 'pymatlab.ScanInput'>
ScanOutput                       type        <class 'pymatlab.ScanOutput'>
filtvec                          function    <function filtvec at 0x1379ff50>
findmovingmotor                  function    <function findmovingmotor at 0x1379fd70>
getSTDMEANfrac_from_startpoint   function    <function getSTDMEANfrac_<...>startpoint at 0x137a01b8>
get_filter                       function    <function get_filter at 0x1379fed8>
get_limits                       function    <function get_limits at 0x137a0050>
get_limits_automatic             function    <function get_limits_automatic at 0x137a0230>
get_limits_channelhist           function    <function get_limits_channelhist at 0x137a02a8>
get_limits_correlation           function    <function get_limits_correlation at 0x137a00c8>
get_limits_corrfrac              function    <function get_limits_corrfrac at 0x137a0140>
h5py                             module      <module 'h5py' from '/reg<...>ython/h5py/__init__.pyc'>
np                               module      <module 'numpy' from '/re<...>thon/numpy/__init__.pyc'>
plt                              module      <module 'matplotlib.pyplo<...>n/matplotlib/pyplot.pyc'>
rdXPPdata                        function    <function rdXPPdata at 0x1379fc80>
runexpNO2fina                    function    <function runexpNO2fina at 0x1379fe60>

In [3]: scaninput = ScanInput()

In [4]: scaninput.fina = "/reg/d/psdm/XPP/xpp23410/hdf5/xpp23410-r0107.h5"

In [5]: scan = rdXPPdata(scaninput)
Reading XPP data from  /reg/d/psdm/XPP/xpp23410/hdf5/xpp23410-r0107.h5
Found pv control object  fs2:ramp_angsft_target
Found scan vector  [ 2800120.  2800240.  2800360.  2800480.  2800600.  2800720.  2800840.
  2800960.  2801080.  2801200.  2801320.  2801440.  2801560.  2801680.
  2801800.  2801920.  2802040.  2802160.  2802280.  2802400.  2802520.
  2802640.  2802760.  2802880.  2803000.  2803120.  2803240.  2803360.
  2803480.  2803600.  2803720.  2803840.  2803960.  2804080.  2804200.
  2804320.  2804440.  2804560.  2804680.  2804800.  2804920.  2805040.
  2805160.  2805280.]
Fetching data to correlate with motor
['IPM1', 'IPM2']
(44, 120, 4)

In [6]: channels = np.concatenate(scan.scandata,axis=0)

In [7]: channels.shape
Out[7]: (5280, 4)

In [8]: get_limits(channels,1,"correlation")
4 channels a 5280 events
indexes that pass filter:  (array([   1,    5,    8, ..., 5266, 5272, 5273]),)
Out[8]:
array([[ 0.00086654,  0.01604564],
       [ 0.67172102,  0.71968567],
       [ 0.00194716,  0.01447819],
       [ 0.80365403,  0.73463468]])

In [9]: plt.draw()

Data visualization with NumPy (arrays) and MatPlotLib (plots).

Saving (and loading) a numpy array (e.g. image) to (from) a file

If you want to save one array (max 2 dimensions), you can use binary numpy file or ascii file:

import numpy as np

# binary file .npy format
np.save("filename.npy", array)
array = np.load("filename.npy")

# txt file
np.savetxt("filename.dat", array)
array = loadtxt("filename.dat")

If you need to save multiple events/shots in the same file you will need to do some tricks (e.g. flatten the array and stack 1d arrays into 2d arrays where axis2 represent event number). Or you could save as an HDF5 file.

You can save an array or several into an HDF5 file (example from pyana):

import h5py

def beginjob(self,evt,env):
    self.ofile = h5py.File("outputfile.hdf5", 'w') # open for writing (overwrites existing file)
    self.shot_counter = 0

def event(self,evt,env)
    # example: store several arrays from one shot in a group labeled with shot (event) number
    self.shot_counter += 1
    group = self.ofile.create_group("Shot%d" % self.shot_counter)

    image1_source = "CxiSc1-0|TM6740-1"
    image2_source = "CxiSc1-0|TM6740-2"

    frame = evt.getFrameValue(image1_source)
    image1 = frame.data()
    frame = evt.getFrameValue(image2_source)
    image2 = frame.data()

    dataset1 = group.create_dataset("%s"%image1_source,data=image1)
    dataset2 = group.create_dataset("%s"%image2_source,data=image2)

def endjob(self,env)
    self.ofile.close()

Or you can group your datasets any other way you find useful, of course.

A comparison with MatLab.

%
%
%
a1 = subplot(121);
loglog(channels(:,1),channels(:,2),'o')
xlabel('CH0')
ylabel('CH1')
a2 = subplot(122);
loglog(channels(:,3),channels(:,4),'o')
xlabel('CH2')
ylabel('CH3')

import matplotlib.pyplot as plt
import numpy as np

a1 = plt.subplot(221)
plt.loglog(channels[:,0],channels[:,1], 'o' )
plt.xlabel('CH0')
plt.ylabel('CH1')
a2 = plt.subplot(222)
plt.loglog(channels[:,2],channels[:,3], 'o' )
plt.xlabel('CH2')
plt.ylabel('CH3')

axes(a1)
hold on
lims(1:2,:) = ginput(2);

axes(a2)
hold on
lims(3:4,:) = ginput(2);

lims = np.zeros((4,2),dtype="float")

plt.axes(a1)
plt.hold(True)
lims[0:2,:] = plt.ginput(2)

plt.axes(a2)
plt.hold(True)
lims[2:4,:] = plt.ginput(2)

fbool1 = (channels(:,1)>min(lims(1:2,1)))&(channels(:,1)<max(lims(1:2,1)))
fbool2 = (channels(:,2)>min(lims(1:2,2)))&(channels(:,2)<max(lims(1:2,2)));
fbool = fbool1&fbool2
loglog(channels(fbool,1),channels(fbool,2),'or')

fbool3 = (channels(:,3)>min(lims(3:4,3)))&(channels(:,3)<max(lims(3:4,3)))
fbool4 = (channels(:,4)>min(lims(3:4,4)))&(channels(:,4)<max(lims(3:4,4)));
fbool = fbool3&fbool4
loglog(channels(fbool,3),channels(fbool,4),'or') 

fbools0 = (channels[:,0]>lims[:,0].min())&(channels[:,0]<lims[:,0].max())
fbools1 = (channels[:,1]>lims[:,1].min())&(channels[:,1]<lims[:,1].max())
fbools = fbools0 & fbools1

fbools2 = (channels[:,2]>lims[:,2].min())&(channels[:,2]<lims[:,2].max())
fbools3 = (channels[:,3]>lims[:,3].min())&(channels[:,3]<lims[:,3].max())
fbools = fbools2&fbools3