pandas on Python kept me from using R

I wanted to analyze this time-series data from GEFCOM visually. The data have missing values and were presented in a way that I didn't want for my analysis.

This looks like a job for R (not numpy/scipy/matplotlib). I had spent a few days playing with R for its supposed excellent visualization and data manipulation capability. I consider myself a Python advocate but Python's main graphics package, matplotlib, is no where near what comes with R. Still, while R's graphics capabilities are great, it wasn't enough for my needs; I needed to work with several visualization packages each with its own way of data input formatting. The all-too-common task of data munging is what I had to deal with. Still, was the data manipulation capability of R enough to keep me on R? After reviewing my analysis needs and what Python had to offer, the answer was: no.

The crucial component in this decision is the existence of the pandas package. In one package, I can read a csv with missing values, reformat it as time series, slice it & dice it any way I want, and plot them for a quick exploratory analysis. Then, in Python syntax, I can reformat the data to feed into another program. When I say, 'Python syntax' I don't mean something like textwritingfunc("some text", anopenfile). I mean a (readable) pseudocode-like syntax that makes use of set and list comprehensions like

out=[ atimeseries:[mathematicaltransformation(index,value) \

for index,value in atimeseries] for atimeseries in alotoftimeseries]

//this could be a python generator to delay processing of the list until file writing to minimize memory usage.

somefile1=openfile("out1.txt")

somefile2=openfile("out2.txt")

//program 1 wants data in a "stacked" format

write header meta data

for atimeseries in mytimeseries: somefile1.writeline(transformedtimeseries in out)

//program 2 wants data in a "wide" format

maybe prog2 doesn't need a header

somefile2.writeline([transformedline for atimeseries in out for alotoftimeseries]) 

This kind of simple yet powerful syntax takes care of the idiosyncratic data input needs of different visualization software. And, if you're a smart coder, you can generalize your code with object-oriented programming for output to multiple programs with minimal fuss. I'm not sure R can do that elegantly.

Conclusions:

• R is a comprehensive statistical tool for analysis and visualization. If all you need to do is in R, stay in R.
• Python is a viable alternative to R for basic to intermediate statistical analysis and visualization. However, expect Python to acquire more statistical tools with time. Even if that's not the case, Python is already well-known as a "glue" language that allows you to pass around data from/to different software thereby scripting an analysis process.
• Python is a better general programming language (it's really hard to say that this is an opinion).

C vs. Python in the Context of Reading a Simple Text Data File

Problem: Using C, read a Matrix Market file generically and versatilely.
Generically: using the same code (at run-time), interpret a line of data in various ways

I'm trying to do this in C. Coming from Python, C makes me feel CRIPPLED. Laundry list of offenses:

• Primitive string operations...oh actually they are just character arrays.
  example: char teststring[10]="asdf"; (OK initialization) teststring="qwer"; (FAIL what I expect is a reassignment fails..turns out you must call string copy and make sure the destination has enough space.).
• No negative indexing
• Can't get length of an array
• Can't handle types as objects. You can't say: if this object is of this type, then do this. Types are 'hardcoded'. You can't return a type.
• Case labels must be known at compile-time. The C compiler can't understand that if I use a function output with a known input at compile time to make a label that it could be known at compile-time. As a mathematically-inclined programmer, all I care about is functionality. I don't like to think about the preprocessor as a separate process.
• Your program can compile and work but in the realm of "undefined behavior". Good luck finding a bug caused by "undefined behavior".

Which brings me to some C items that have overlapping functionality: enum, arrays, and X-Macros. They can all be viewed as lists, but macros live in the preprocessor; you cannot iterate over enums and they are essentially only a list of integers; and if you want to know the size of an array, you can only do it at compile-time with a macro hack.  So if you need functionality on the same list that had to be implemented in the preprocessor and at run-time, you have to violate DRY.

While I do think lower-level functions are essential for some purposes, this exercise makes me appreciate Python more (yes particularly Python and not other interpreted languages). In Python, I can read the lines in as a matrix of strings, then I could vectorize type conversion operations in a few lines of code.

So, I've resorted to code generation to create functions for all cases which solves half the problem. I won't bother trying to make things more generic. If it could be done, I expect that it's going to look messy.

4/17 Update: I've finished my reader. The "main()" file is here. I think I organized my procedure well but it took alot of code. It's ridiculous. I could accomplish the same with acceptable speed in about two dozen lines of Python code. I know this because I've written code to read data generically in Python several times.

Mathematica vs. Maple vs. SAGE/scipy UX

This post comes from significant experience with Maple and numpy/scipy/python.

I've been using Maple and numpy for a few years. I stopped upgrading my Maple at version 13 because improvements didn't compel me to pay for a newer version for my use. I chose Maple as a happy medium between the numerical and matrix-focused MATLAB and the symbolically-focused Mathematica. They can all do symbolics and matrix operations but they differ in how elegantly they're performed. But, as I gave Maple more complex mathematical tasks such as the calculus involved with Fourier transforms, and symbolic tensor & matrix operations, it became clear that it wasn't up to the task. I knew Mathematica was superior for those tasks but I chugged along with Maple until I got to George Mason where they had a site license for Mathematica.

Right off the bat I could see a clear superiority in the user experience (mind that I didn't use Maple beyond ver 13). Vectors and matrices are just lists. In Maple there are objects called vectors, lists, tables, and matrices and you sometimes had to import them from a module. Clearly, in just this area, Mathematica is superior. Mathematica includes alot of functionality out of the box. In contrast, in Maple, I had to import modules for what I considered basic functionality which makes the experience not as consistent.

In comparison to other math environments:

• SAGE: SAGE uses a web interface which can't come close to a desktop application's functionality. The UI does not come close to Mathematica. Also, the open-source symbolic packages it uses are primitive compared to Maple let alone Mathematica. Oh and good luck installing it on Windows.
• numpy/scipy: Numpy and scipy give mathematical functionality as opposed to being a math environment; basic ones at that. But, its indexing coolness does exist in Mathematica.
• Mathematica has pattern-matching..others don't enough said.

For math-centric programming, Mathematica is a high bar to get to in terms of consistency, documentation, functionality, and general experience. My general gripe about open-source is the lack of an umbrella vision that moves a project in a certain direction. SAGE is the best (the only) hope but it has a long way to go. Even in the best scenario, getting different python packages seamlessly working together is doubtful (eg. look at my sym2num function). Conclusion: If all you need is basic functionality, Octave, SAGE, scipy, sympy, and maxima could suit you. For more advanced tasks, be prepared to pay up.

PS: This post will change as I learn more.
Support for my opinion.

Personal Note: Adjustment in Career Path Towards Data Science

I want this blog to be about the work that I do and not about my person. But there will be a change in the the content of the blog which can only be explained by the change in my professional situation. So far, computational topics have been about FDTD simulation, molecular dynamics simulation, python, and scientific computing in general. I've been a graduate student at Vanderbilt University for four and a half years where I've learned an incredible amount in the areas of nanoscale science & engineering and computing, both of which are areas I did not have the background for. This was all under a mechanical engineering degree but I'm leaving with 'just' a master's degree in "mechanical engineering" as I did not pass my PhD qualifying exam in...(classical) mechanical engineering. Now, I was quite disappointed for a while since I enjoyed my research very much and was looking forward to making contributions to energy conversion device research using computation.

But, looking back only a few months later, I'm glad it happened. I felt that I had a programming talent lurking inside me as I was conducting my computational research and it would be wise to enhance my computational skills which would qualify me very well for certain job types. So, I've started a (second) master's in computational science at George Mason University customizing my curriculum for 'big data'. I've identified 'big data' as a field that I want to get into because I feel that I have the right combination of communication, business, and technical skills to be successful in this field. This potentially means that I might turn my back to science and engineering as an application area of my skills which is something a bit discomforting to me as I've always seen myself as a scientist and engineer. At the same time, advanced scientists and engineers have the potential to be great data scientists.

What does this mean for the blog? I still expect to post about (pure) scientific computing but not about applications like finite-difference time-domain simulations and molecular dynamics simulations. Python will stay with me as a base of my programming skill but it will be complemented by other computing languages that I will learn like SAS, R, C, and Fortran (working on working in HADOOP!). This is what I'm able to predict for now.

Ending with a positive personal note, for most of the past decade I've been striving to get into a job type that I want. I'm very optimistic that I'm almost there.

Making My Near-Field Radiation Calculator Application for nanohub.org

I wanted to share a bit of my experience developing my python application for nanohub.org. For me, there were two main issues: One is configuring the nanohub computing environment to resemble my development environment (and you'll need to figure out the details of installing your tool by looking at other tools and some docs that can help). The other is using Rappture to implement the functionality I want for my application.

To set up my development environment in nanohub, I first delved into how I could make a self-contained python installation. This had me looking into pythonbrew and virtualenv with pythonbrew being more useful. But the admins at nanohub were willing to install the packages that I needed so I didn't have to mess with that too much. Nonetheless, these two projects are great for creating portable python projects.

After I took care of the development environment, I went ahead and started learning about Rappture. I needed some flexibility and what I call 'dynamic input'. Now I could have made a simple application but then it wouldn't be as useful. I had to bend the API to suit my needs. I ended up merging what should have been two simple applications into one more complicated application. That is: I had two sets of inputs for the same calculation. I'll list the issues I had with Rappture:

• Lacking documentation especially for the Python API. This is compensated for by looking at the API for the other languages and the many examples. I mainly figured out how Rappture works by example.
• Lack of computation status output as addressed by this page.
• Not all GUI elements are in the GUI GUI builder.
• Could not generate an authentic contour plot. I used the elements field and cloud to generate a contour plot. This probably wasn't the intended use of the of the elements but it worked...horribly. It struggled to display just a 50x50 grid calculation not to mention that I couldn't put axis labels. Also of general note, elements take strings to input numbers which is inefficient (but works). In my case, I wanted to be able to display at least 1000x1000 points. My solution was to ditch this GUI element and just use an image generated by matplotlib's contourf.
• 'Static' input. I wanted a user of my program to be able to change and check inputs in certain ways.
• I wanted the user to be able to add a selection to a drop-down list. Could not be done. Instead, I have the user reviewing a table of named text input so that the name can go into a text box which should have been a selection.
• Have an input value checked based on another input value. Can be done functionally by writing code that gets processed after hitting the 'Simulate' button. I don't know how I could hook into the GUI input elements so that they can be checked before the inputs are processed for simulation. For example, in my program I wanted to check that the value of 1.0 wasn't between two input values.
• Limited selection of units. (eg. no rad/s.) There should be functionality so that any unit can be input and can be converted to a compatible unit.

Despite these concerns, I believe Rappture is great for those simple (from an input/output point of view) programs that were originally designed for textual input and output. For more complicated programs, Rappture poses limitations. But then there are other developed and more sophisticated GUI APIs out there. So a better Rappture would be duplicating efforts of those other APIs.

Here's how it works:
First you get a page where you can enter named electric permittivity functions.

Then you input the simulation parameters and options.

Then you get various outputs. I'm showing the cool ones here:

The characteristic nearly monochromatic near-field radiation between SiC and SiC

The above graph is an integration over the y-axis of the following function in two variables.

Task 6: Basline Near-Field Calculation

I'm going to be doing FDTD stuff but need a baseline transport calculation using a closed-form expression. Maple, Maxima, Sympy couldn't do the integration. So I had to do it numerically in a more manual way. So I used my trusty Python environment with the following modules that did most of the work:

- sympy: I used this package to build the final complicated expression that gives the heat transfer. Unfortuately, it's immature software. For example it couldn't solve: sympy.solve((x**2)**0.5-9). Also the lambda features were confusing since there were two lambda features in the package.
- numexpr: 1.4 I used this blazingly fast expression evaluation package to evaluate the built expression. (I had to comment line 263 in necompiler.py because I was working with imaginary numbers. It's an optimization step that needs to order numbers and there is no order to imaginary numbers) See better workaround.

11/8: It's published on nanohub.org.

Prototype GUI built using hubzero's Rappture GUI builder

The code for the core of the program is below.

#Author: Majid S. al-Dosari

#http://www.thermalfluidscentral.org/encyclopedia/index.php
#/Near-field_radiative_transfer_between_two_semi-infinite_media
#Calculate Bulk 1-D Near-Field Radiation between two Bulks.
#Input:
#1. permittivity of both bulks
#2. distance b/w them
#3. Temperature of each bulk T1,T2
#Output:
#Power as a function of frequency or just total power

#physics notes
#the integrand is split into TE and TM components
#these components are further split into propagating (omega/c<1 data-blogger-escaped-and="and" data-blogger-escaped-c="c" data-blogger-escaped-evanescent="evanescent" data-blogger-escaped-omega="omega" data-blogger-escaped-waves="waves">1)

#computation procedure
#1. build expression for TE(S-polarization) and TM(P-polarization) waves.
#"s" in the web reference.
#See 'sept' and 'sest' names below as examples
#2. convert sympy expression to numexpr
#3. use evalq for total power, evalqw for power as a function of frequency,
#evalqwsections to get power as a function of frequency when the practical
#upper limit of integration in wavevector is unknown (infinity)
#they share the common argument sequence
#(stringformofintegrand,beta1,beta2,T1,T2,nb=1000,omega1=0,no=1000)
#beta here is normalized to omega/c. omega2 is determined as that frequency at 
#which 99.99...?% blackbody radiation integral is covered at the larger 
#of T1, and T2.
#nb and no are the number of points discretized on beta and omega respectively
#don't include beta=1 in the range b/w beta1 and beta2.
#use beta2=.9999 as an upperlimit for integrating propagating waves and
#use beta1=1.0001 as a lower limit for integrating evanescent waves.
#for evalq and evalqwsections, if beta2 is in the evanescent range (>1), it will
#be taken as an integration 'chunck' size and will continue integration until
#asymptotes. be careful, because if you know TM waves's contribution 
#are negligible it may continue to integrate to ridiculous wavevectors and 
#therefore give unrealistic results. there is an absolute limit for beta2
# around a few interatomic distances
#programdata is a dictionary that holds frequencies and wavevectors that the
#integrand was evaluated over. i use them for plot axes.

#example: evalq(sym2num(str(sept)),1.001,50,300,299,no=3e3,nb=3e2)

#use the output from evals to make a pretty picture of the evaluated func
#eg use matplotlib's contourf(swapaxis(evals output))


#coding note:
#a=1+2j
#numexpr.evaluate('(a*1j)+0')
#gives TypeError: no ordering relation is defined for complex numbers
#commented out line 263 #ordered_constants.sort() in necompiler
#see 'HACK' in sym2num

#evaluate new option to sub output inplace of input

import numpy
import sympy
import numexpr


#might be able to get away w/ 100sx100s grid
#.0001 and .9999 , 1.00001


#Material Permittivity Input
#---
def invcm2rad(ic):return ic*2*numpy.pi*2.998e10
consts={'c':2.998e8,'hbar':1.054571596e-34,'kb':1.380650277e-23
,'omega_LO':invcm2rad(969),'omega_TO':invcm2rad(793),'Gamma':invcm2rad(4.76)
,'eps_inf':6.7}
eps_inf,omega_LO ,omega_TO ,Gamma,omega=\
sympy.symbols('eps_inf omega_LO omega_TO Gamma omega')
eps_SiC=eps_inf*(1+(omega_LO**2-omega_TO**2)\
        /(omega_TO**2-1j*Gamma*omega-omega**2))

#eps_Au=1-(13.71e15**2)/(omega*(omega+2*numpy.pi*1j*4.05e13))

#indexed epsilon
eps=[1,eps_SiC,eps_SiC]
#eps=[1,eps1,eps2]



beta,c=sympy.symbols('beta c')
beta=beta*(omega/c) #in terms of omega/c
def gama(j):return (eps[j]*((omega/c)**2)-beta**2)**.5
#def gama(j):return (eps[j]-(beta/(omega/c))**2 )**.5 #in terms of omega/c
#0j to fix ValueError: negative number cannot be raised to a fractional power
#not needed with sympy i think

def r0_s(j):return (gama(0)-gama(j))*(gama(0)+gama(j))
def r0_p(j):return (eps[j]*gama(0)-gama(j))/(eps[j]*gama(0)+gama(j))
r0=dict({'s':r0_s,'p':r0_p})

d=sympy.Symbol('d')
def sevan(sorp):
    return sympy.im(r0[sorp](1))*sympy.im(r0[sorp](2)) \
        *beta*sympy.exp(-2*sympy.im(gama(0))*d)/ \
    sympy.abs(1-r0[sorp](1)*r0[sorp](2)*sympy.exp(-2*sympy.im(gama(0))*d))**2
#the sympy on nanohub doesn't have sympy.Abs rather it has sympy.abs
def sprop(sorp):
    return (1.0/4)*beta*(1-sympy.abs(r0[sorp](1))**2)*(1-sympy.abs(r0[sorp](2))**2)\
    /sympy.abs(1-r0[sorp](1)*r0[sorp](2)*sympy.exp((2j)*gama(0)*d))**2
#the s funcs /give/ the expr. they aren't funcs themselves

hbar,T,kb=sympy.symbols('hbar T kb')
Thetaexpr=hbar*omega/(sympy.exp(hbar*omega/(kb*T))-1)
Thetaexpr=Thetaexpr.subs(consts)
def Theta(omega,T):
    if T==0: return 0
    return Thetaexpr.subs({'T':T,'omega':omega})

def q_omega(omega,T1,T2,intofs):
    return ((1./sympy.pi)**2)*(Theta(omega,T1)-Theta(omega,T2))*intofs

#import scipy
#from scipy import integrate
#useless
def intofs(sfunc,beta1,beta2):
    #sfunc1=lambda b: sfunc(b+0j)#needed for numexpr
    #return scipy.integrate.quad(sfunc1,beta1,beta2)[0] #slowwwww
    #i made this func bc i could use scipy quad or something else
    return sympy.mpmath.quad(sfunc,[beta1,beta2])



#compress const array

#evaluated numerically from the start in a simple way
def nintofs(numstr,beta1,beta2,omega,num=100):#sfunc w/o omega :(
#when i subs in an omega i get weird results
    [beta,omega],db,do=geninputarray(beta1,beta2,num,omega,1234,1)
    return db*numpy.sum(numexpr.evaluate(numstr))
    

def evals_(numstr,*args,**kwargs):#sfunc w/o omega :(
#when i subs in an omega i get weird results
    [beta,omega],db,do=geninputarray(*args,**kwargs)
    return (numexpr.evaluate(numstr)),db,do #with zeros you get NaNs
    #swap(reshape())
def evals(sstr,*args,**kwargs):
    so,db,do=evals_(sstr,*args,**kwargs)
    so=numpy.reshape(so,(args[2],args[5]))
    #swap args for reshape other way for contourf, then swapaxis
    return so,db,do
    

def evalDTs(DTstr,omegas):#sfunc w/o omega :(
#when i subs in an omega i get weird results
    omega=omegas
    if omega[0]<.1:omega[0]=.1 #get's rid of NaNs
    return numexpr.evaluate(DTstr)

#the limits of integration
planckexpr=(hbar*omega**3/(4*sympy.pi**3*c**2))\
    /(sympy.exp(hbar*omega/(kb*T))-1)
planckexpr=planckexpr.subs(consts)
def planckd(o,T):return planckexpr.subs({'omega':o,'T':T})
#dplanckexpr=planckexpr.diff(omega)

#from scipy import optimize
def omegalim(T): return consts['c']*2*numpy.pi/(1000e-6/T) #-.000321+1 of BB
#% this might have an effect

def solvebeta(omega): return omega/consts['c']
    #return sympy.solve(gama(0).subs(consts).subs({'omega':omega}))[0]
    #it can't solve this!

def betalim(T): return omegalim(T)/consts['c']

def evalqwint(sfunc,beta1,beta2,T1,T2,omega1=1e-20,no=500):
    os=numpy.linspace(omega1,omegalim(max(T1,T2)),no)
    sfb=lambda o:sfunc.subs({'omega':o})#gives func of sympy
    spec=[intofs( lambda b:sfb(o).evalf(subs={'beta':b})  ,beta1,beta2) for o in os]
    #spec=[nintofs( sym2num(str(sfunc)) ,beta1,beta2,o,num=1e3) for o in os]
    #np.random.seed(123)    
    #DTs=np.random.rand(len(os))
    #this is introducing a good amoutn of error in SiC for d=100e-9
    DTs=(evalDTs(sym2num(str((Theta(omega,T1)-Theta(omega,T2)))) ,os))
    qw= spec*DTs*(os/consts['c'])/numpy.pi**2# no 2pi ??
    return qw,os[1]-os[0] #dw


def evalqw(sstr,beta1,beta2,T1,T2,nb=1000,omega1=0,no=1000):
    #evals(sym2num(str(sppt+spst)),0,.999999,1000,0,omegalim(300),1000)
    so=evals(sstr,beta1,beta2,nb,omega1,omegalim(max(T1,T2)),no)
    #so=evals(sstr,beta1,beta2,nb,1.7e14,1.85e14,no)
    os=numpy.real(programdata['os'])
    dBetas=so[1]*(os/consts['c'])
    #np.random.seed(123)    
    #DTs=np.random.rand(len(os))
    #this is introducing a good amoutn of error in SiC for d=100e-9
    DTs=(evalDTs(sym2num(str((Theta(omega,T1)-Theta(omega,T2)))) ,os))
    qw= numpy.sum(numpy.real(so[0]),axis=0)*DTs*dBetas/numpy.pi**2# no /2pi ??
    return qw,so[2] #dw


import scipy
from scipy import integrate
def evalq(*args,**kwargs):
    if args[2]<=1: o=evalqw(*args,**kwargs) #for propagating waves
    else: o=evalqwsections(*args,**kwargs)
    #return sum(o[0])*o[1]
    return scipy.integrate.trapz(o[0],dx=o[1])


def evalqwsections(*args,**kwargs):
    #take beta2 as sectioning
    #if beta2>1:
    args=list(args)
    interval=args[2]-args[1]
    ac=asymptote()
    passn=1
    xlim=(numpy.pi/.5e-9)/(omegalim(max(args[3],args[4]))/consts['c'])#~=7e3
    #(2*3.14/10e-10)/ #per zhang review of NF thermal rad    
    #say 10angstrom is local limit
    while 1:
        if args[1]>xlim:break
        o=evalqw(*args,**kwargs)
        if passn==1:cur=numpy.zeros(len(programdata['os']))
        cur+=o[0]
        sm=numpy.sum(o[0])*o[1]
        if ac.didit(sm)==True:break
        passn+=1
        #shift for next interval
        args[1]+=interval;args[2]+=interval;
    return cur, o[1]


class asymptote():
    #expectation of the trend is that it goes up then goes down in the end
    #it could vary up and down in b/w
    def __init__(self):
        self.nums=[]
        self.curmax=float("-inf")
    def didit(self,num):
        tol=.01
        if len(self.nums)==0:self.curmax=num
        self.nums.append(num)
        if self.up()==True:
            self.curmax=num
            return False
        if self.up()==None: return None
        if self.up()==False and self.curmax*tol>num:return True
        else: return False
    def up(self):
        if len(self.nums)<2:return data-blogger-escaped-if="if" data-blogger-escaped-none="none" data-blogger-escaped-self.nums="self.nums">0:return True
        else: return False

onejay=1j
def sym2num(astr):
    """
     to go from a sympy expr to a numexpr expr
    """
    astr=astr.replace('I','1j')
    astr=astr.replace('1j','onejay') #A HACK!!!
    astr=astr.replace('im','imag')
    astr=astr.replace('re','real')
    astr=astr.replace('Abs','abs')
    return astr
 

#a hold for intermediate data
programdata={}

def geninputarray(betamin,betamax,bnum,omegamin,omegamax,onum):#,num=None):
#onum=100,bnum=100,num=None
#    ,betamin=1e-20,betamax=2e6,omegamin=1e-20,omegamax=2e14):
    """inputs in numexpr
    extreme:say you want just one line instead of 2D:
        specify a close bmin and bmax with bnum=1
        betamin will always be evaluated. to evaluate at a certain beta put it
        in betamin
    """
#    if num!=None:#if num specifided
#        bnum=onum=num
    #bs=numpy.add(numpy.linspace(betamin,betamax,num=num),0j)
    bs=numpy.zeros(bnum,dtype='c8') #could use bigger float but i dont care
    bs+=numpy.linspace((betamin),(betamax),num=bnum,endpoint=False)
    #add 0j needed to work in expr
    os=numpy.zeros(onum,dtype='c8')
    os+=numpy.linspace((omegamin),(omegamax),num=onum,endpoint=False)
    if os[0]==0:os[0]=1e-20
    if bs[0]==0:bs[0]=1e-20
    programdata.update({'os':os,'bs':bs})
    return (cartesian([bs,os]).transpose()
        ,(float(betamax)-betamin)/len(bs),(float(omegamax)-omegamin)/len(os))
        
    #fxy eval, dx*dy

#evanescent waves test expressions
sep=sevan('p').subs(consts)
sept=sep.subs({'d':100e-9})
ses=  sevan('s').subs(consts)
sest=ses.subs({'d':100e-9})
#propagatin waves test expressions
spp=sprop('p').subs(consts)
sppt= spp.subs({'d':1})
sps=sprop('s').subs(consts)
spst=sps.subs({'d':1})



import numpy as np
#WOW!
def cartesian(arrays, out=None):
#http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays
    """
    Generate a cartesian product of input arrays.

    Parameters
    ----------
    arrays : list of array-like
        1-D arrays to form the cartesian product of.
    out : ndarray
        Array to place the cartesian product in.

    Returns
    -------
    out : ndarray
        2-D array of shape (M, len(arrays)) containing cartesian products
        formed of input arrays.

    Examples
    --------
    >>> cartesian(([1, 2, 3], [4, 5], [6, 7]))
    array([[1, 4, 6],
           [1, 4, 7],
           [1, 5, 6],
           [1, 5, 7],
           [2, 4, 6],
           [2, 4, 7],
           [2, 5, 6],
           [2, 5, 7],
           [3, 4, 6],
           [3, 4, 7],
           [3, 5, 6],
           [3, 5, 7]])

    """

    arrays = [np.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = np.prod([x.size for x in arrays])
    if out is None:
        out = np.zeros([n, len(arrays)], dtype=dtype)

    m = n / arrays[0].size
    out[:,0] = np.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian(arrays[1:], out=out[0:m,1:])
        for j in xrange(1, arrays[0].size):
            out[j*m:(j+1)*m,1:] = out[0:m,1:]
    return out