A beginning Python programming challenge: replace slow array.index() calls.

For simplicity, consider the following arrays:

i      0  1   2  3  4  5  6  7  8  9  10

order1 6  9   7 10  5  8 3  2  4   0   1

order2 7  6  10  8  5  9 1  4  2   0   3

Now consider the following Python snippet of code:

myindex = range(11)

for i in range(11):

    myindex[i] = order2[order1.index[i]]

As an example consider a current i value of 0. The value of index1[0] is 6. The position of 6 in order2 is of 6 in index2 array is 1(first element has index0). Therefore myindex[0] = 1.

On the other hand, the value of index1[1] is 9. The ordinal position of 9 in index2 is 5. Hence myindex[1] = 5.

Problem 1. Can you replace the code without using dictionaries and the index function?
You may not modify inplace the order1 and order2 arrays. You may create an additional extra array.

Problem 2.  Rewrite the code to use a dictionary. This is much easier.

Imagine your are doing the for loop of 1,000,000 array items! Have fun!

Python: Generating all subsets of size k for n objects.

In our main blog, we presented a routine for generating all subsets of set, expressed as an array or vector of indices to the n objects (0 to n-1).See Generating all subsets of a set with n objects In this article we present an "odometer" method of generating all subsets of specific size k where k < n. For example, for n = 5 and k = 3, we expect the function to return the set containing the following vectors: [0,1,2], [0, 1,3], [0,1,4], [0,2,3] and so on up to [2,3,4]. Here is our Python code.

# -*- coding: utf-8 -*-



"""

File      ksubsets.py

Author    Ernesto P. Adorio, Ph.D.

UPDEPP, UP Clarkfield  

Universty of the Philippines Extension Program in Pampanga

Desc      generates subsets (as indices) with k elements from a set of size n.

Version   0.0.1  2011.07.03

"""





def  ksubsets (n,  k):

"""

generates the list of subsets (as index array) of n elements

chosen k at a time. Combinations actually!

"""

#sanity check

if n <0 or k <0 or n - k < 0:
        return None
    N = range(k)
    retval = [N[:]]
    print "initial N=", N
    
    #update current N vector.
    level = k-1
    while level >=0:

d = N[level]

#test if element is already  at maximum value.

if d < n - k + level:  
            # No, increment the value.  
            N[level]+= 1
            #repair trailing elements.
            for i in range(level+1,  k):
                N[i] = N[i-1] + 1
            #append newly generated subset.
            retval.append(N[:])
            level = k-1
        else:
            # position N[level] already at maximum value!
            level -= 1
    return retval
    
for i,  subset in enumerate(ksubsets(5,  3)):
    print i,  subset

When the above program is run, it outputs the expected results:

0 [0, 1, 2]
1 [0, 1, 3]
2 [0, 1, 4]
3 [0, 2, 3]
4 [0, 2, 4]
5 [0, 3, 4]
6 [1, 2, 3]
7 [1, 2, 4]
8 [1, 3, 4]
9 [2, 3, 4]

How does the code work? First it initializes the first member of the set of subsets as
[0, 1, ..., k-1]. Next, starting from the rightmost (level = k-1)position it checks whether the current element is at maximum. If it is not it is incremented and any trailing positions adjusted to be 1 greater than at the preceding position. Otherwise the level is decremented and the process repeats.

Exercise: write a generator version of the above code. We will present the solution on the next revision.

The itertools module has a combinations function. After importing, one simply writes:
list(combinations (range(5), 3)).. Python just makes it TOO easy!

Python in Electrical Engineering: Converting between Y and Delta impedances

Given the branch impedances in delta configuration, what is the equivalent circuit impedance if connected in a Y configuration? the folowing short Python codes ease the conversion between Y and Delta impedances.

"""
file       deltay.py
author     Ernesto P. Adorio, Ph.D.
           UPDEPP (U.P. at Clarkfield)
desc       translate delta impedances to Y impedances.
version    0.0.1    04.21.2011
"""



def delta2yZ(Zab ,Zbc, Zca):
    """
    translate delta impedances to Y impedances.
    Arguments
       Zab, Zbc, Zca - complex impedances of Deta connection.
    Output
       Za, Zb, Zc - equivalent Y connection impedances to commom point.
    """
    denom = Zab + Zbc + Zca
    return (Zca * Zab/denom,Zab* Zbc/denom, Zbc * Zca/ denom)

def y2deltaZ(Za, Zb, Zc):
    """
    translate Y impedances to delta impedances.
    Arguments
       Za, Zb, Zc - complex Y connected impedances 
    Output
       Zab, Zbc, Zca - equivalent delta connection branch impedances.
    """
    num = Za*Zb + Zb*Zc + Zc* Za
    return (num/ Zc, num / Za, num/Zb)

We will continue this with an example later.

Exercise: State the formulas used by the above routines y2deltaZ(), delta2yZ() routines. It is so easy to write the functions that we did not bother to write the documentation for this first version, but we will later!

Python in Elecrical Engineering: Computing Symmetrical Components.

Given three-phase voltages V1, V2, V3 in phasor (r, phi) form, where r is the magnitude of the phasor, and phi is the phase angle. the following formulas compute the corresponding symmetrical components:

Null sequence component: $$E_0 = (1/3)(V1 + V2 + V3)$$.

Positive sequence component: $$E_1 = (1/3)(V1 + aV2 + a^2 V3)$$.

Negative sequence component: $$E_2 = (1/3)(V1 + a^2V2 + a V3)$$.

Here the a is the operator which rotates the phasor by 120 degrees. Power engineers usually specify the input voltages or current in phasor form. The formula shows we have to add phasors termwise and thus the routine to compute symmetrical components first convert them to complex numbers. Complex numbers are already handled naturally in Python! and here is the full source code:

"""
file    symcomp.py
author  Ernesto P. Adorio, Ph.D.
UPDEPP (U.P. Clarkfield)
desc    basic symmetrical components for three phase phasors.
Phasors are tuples of the form (r, phi) where r and phi are the
magnitude(absolute value) and phi is the phase angle.
version 0.0.1 april 20, 2011
"""


from math  import *
from cmath import *

zTOL = 1.0e-8
DTOR = pi / 180.0
RTOD = 180.0 / pi



ONETWENTYRAD = 120 * DTOR
_a_  = rect(1, 120 * pi/ 180.0)
_a2_ = _a_.conjugate()

def dtor(degree):
"""
Converts degree to radians.
"""
return degree * DTOR

def rtod(radian):
"""
Converts radians to degree.
"""
return radian * RTOD

def  z2pair(z):
"""
returns z in pair form.
""" 
return (z.real, z.imag)

def  pair2z(re, im):
"""
returns a complex number from the components. 
"""
return re + im*1j

def  polar2pair(v):
"""
(r, phi) to (re, im) pair form.
"""
# Special checking for zero imaginaries!
# some computations result in an additional 0j
r, phi = v
if type(r) == type(1j):
r =r.real
if type(phi) == type(1j):
phi = phi.real

z = rect(r, phi)
return z.real, z.imag

def  pair2polar(re, im):
return polar(pair2z(re,im))

def  zround(v):
if abs(v[0]) < zTOL:
        return (0.0, 0.0)

def  a(v):
     """
     Applies the a operator to a phasor v in polar form.
     It adds a 120 degree to the phase of the phasor v.
     if v is real, the angle is zero.
     """
     if type(v) != type((0,0)):
        v = (v, 0)
     r     = abs(v[0])
     theta = v[1]

     newangle = theta + ONETWENTYRAD
     return (r * cos(newangle), r * sin(newangle))

def a2(v):
     """
     Applies the a operator to a phasor v in polar form.
     It adds a 120 degree to the phase of the phasor v.
     """
     if type(v) != type((0,0)):
        v = (v, 0)
     r     = abs(v[0])
     theta = v[1]

     newangle = theta + 2*ONETWENTYRAD
     return (r * cos(newangle), r * sin(newangle))


def  symcomp(v1, v2, v3):
     """
     Returns the symmetrical components of the three phasors v1, v2, v3
     which are in tuple (r, theta) form.
     """
     #Convert first to complex rectangular form.
     v1z  = rect(v1[0], v1[1])
     v2z  = rect(v2[0], v2[1])
     v3z  = rect(v3[0], v3[1]) 

     av2  = _a_ * v2z
     a2v2 = _a2_ * v2z
     av3  = _a_* v3z
     a2v3 = _a2_* v3z

     #Null sequence  component.
     E0 = polar((v1z.real+ v2z.real+  v3z.real)/3.0 +  (v1z.imag+ v2z.imag+ v3z.imag)/3.0*1j)

     #Positive sequence component.
     E1 = polar((v1z.real + av2.real +  a2v3.real)/3.0+ (v1z.imag+ av2.imag+ a2v3.imag)/3.0*1j)

     #Negative sequence component.
     E2 = polar((v1z.real + a2v2.real +  av3.real)/3.0+ (v1z.imag+ a2v2.imag+ av3.imag)/3.0*1j)

     return (E0, E1, E2)


def  symcompz(v1, v2, v3):
     """
     Returns the symmetrical components of the three phasors v1, v2, v3
     which are in tuple (r, theta) form.
     """
     av2  = _a_ * v2
     a2v2 = _a2_ * v2
     av3  = _a_* v3
     a2v3 = _a2_* v3

     #Null sequence  component.
     E0 = polar((v1.real+ v2.real+  v3.real)/3.0 +  (v1.imag+ v2.imag+ v3.imag)/3.0*1j)

     #Positive sequence component.
     E1 = polar((v1.real + av2.real +  a2v3.real)/3.0+ (v1.imag+ av2.imag+ a2v3.imag)/3.0*1j)

     #Negative sequence component.
     E2 = polar((v1.real + a2v2.real +  av3.real)/3.0+ (v1.imag+ a2v2.imag+ av3.imag)/3.0*1j)

     return (E0, E1, E2)

def  symcomp2phasors(E0, E1, E2):
     """
     Recreates the phasors form the symmetrical components.
     """
     V1 = polar(rect(E0[0], E0[1]) + rect(E1[0], E1[1]) + rect(E2[0], E2[1]))
     V2 = polar(rect(E0[0], E0[1]) + _a2_* rect(E1[0], E1[1]) + _a_ *rect(E2[0], E2[1]))
     V3 = polar(rect(E0[0], E0[1]) + _a_* rect(E1[0], E1[1]) + _a2_ * rect(E2[0], E2[1]))
     return V1, V2, V3


if __name__ == "__main__":
   #extreme cases.
   I1 = polar(10)
   I2 = polar(0)
   I3 = polar(0)
   print symcomp(I1, I2, I3)   
  
  

   #extreme cases, balanced system.
   i1 = 1
   i2 = -0.5+sqrt(3)/2.0 * 1j
   i3 = -0.5-sqrt(3)/2.0 * 1j

   I1 = polar(i1) 
   I2 = polar(i2)
   I3 = polar(i3)
   E0, E1, E2 = symcomp(I1, I2, I3)   
   print "original phasors=", I1, I2, I3
   print "symmetrical components:", E0, E1, E2
   phasors = symcomp2phasors(E0, E1, E2)
   print "recovered phasors:", phasors
   #include more here! from published books or other sources.

We have not yet fully tested the above code, and we would appreciate it if our readers point out any any errors in the program. When the current version 0.0.1 is run, it outputs the following:

$ python symcomp.py
((3.3333333333333335, 0.0), (3.3333333333333335, 0.0), (3.3333333333333335, 0.0))
original phasors= (1.0, 0.0) (0.99999999999999989, 2.0943951023931957) (0.99999999999999989, -2.0943951023931957)
symmetrical components: (7.4014868308343765e-17, 3.1415926535897931) (1.295260195396016e-16, 0.0) (0.99999999999999989, 0.0)
recovered phasors: ((1.0, 9.0639078007724287e-33), (0.99999999999999978, 2.0943951023931953), (0.99999999999999978, -2.0943951023931953))
$

Note that in the second example, the null-sequnce component is the zero vector (zero magnitude). Out routines are able to recover the original phasors from the computed symmetrical componsnts.

Using javascript for data visualization and animations.

Visit neurosoftware.ro.


We are still learning how to run jsxscript in our Wordpress blog.

Adding commas to integer strings for readability

There ase many ways to make more human readable the string representation of an integer.
Click on Activestate.com thousands separator

As an example, 1234567890 looks more humanly readable as 1,234,567,890.

We normally use other people's tested and understood code to save time, but this is one good programming problem.

The manual way to do this is to insert commas before blocks of contiguous three digits.
Strings are supposed to be immutable in Python by design for simplicity. Every changes to a string must be done by a copy. So we convert a string to a list just so to insert characters!

def  withcommas(s):
     p = len(s) -3
     s = list(s)
     while p > 0:
        s.insert(p, ",")
        p -= 3
     return "".join(s)

But suppose that s starts with a '+" or "-" sign? Then the leftmost position for a comma is at the second character or index 1. We add a conditional test so that the code will not be put off by the leading sign.

def  withcommas(s):
     p = len(s) -3
     s = list(s)
     if s.startswith("+") or s.startswith("-"):
        firstpos = 1
     while p > firstpos:
        s.insert(p, ",")
        p -= 3
     return "".join(s)

and finally we include a trivial programming problem for the reader. Instead of 3, let the caller of withcommas() specify a blocklength to the function. We are just getting ready for the time when we will concentrate on creating mathematical tables.

Our aim is to create simple, working first prototypes. The reader again is encouraged to read and search the cited reference above for faster solutions.

Pickled objects in Python

The inventor of Python has given a catchy name to saving objects in files as "pickled objects".
I often use the cPickle module and one can get by using only two powerful functions load and dump in applications.

1. Save a named object to a pickle file.

import cPickle
a = "The quick brown fox"
b = "Hello world."
c = 123

D = { "a": a, "b": b}

Note that we create a dictionary containing the assigned names of the objects, and their corresponding values. An object is stored in external file and it is possible to recover their contents (seconds, minutes, hours, days, even years later)! Now here is the code for storing in the pickfile.

pickfilename = "tmp.pck"
pickfile = open(pickfilename,"w")

cPickle.dump(D, pickfile)


The cPickle.dump() function expects an object and an open file in write mode. If you dump the contents of tmp.pck you will get something like this displayed on the screen:


The idea of using a dictionary is the most natural way of saving a set of variables you are interested at to save and look for values later. Of course if you only wanted to save the variable a, you can do cPickle.dump(a, pickfile)

The name of all active local variables can be accessed by the function vars() which is a dictionary. But it also includes loaded modules which you may not want to be included!

2. Load saved object in pickle file

Nowt that you have created a pickfile, "tmp.pck", let us take a look at it. Dumping it on the console gives

(dp1
S'a'
S'The quick brown fox'
p2
sS'b'
S'Hello world.'
p3

So at least we are at ease that the values are indeed stored for later access.

myvars = cPickel.load(open(pickilename,"r"))

We knew that the name of the dictionary variable stored was "D". But you dont have to use D.
Instead above, we use a new label, myvars.

Now we can recover the values of the former variables by indexing using 'a' and 'b' as keys.Thus,


a = myvars["a"]
b = myvars["b"]


Notes.

The cPickle functions has an additional parameter calle 'protocol'. At this point in time, we refer the interested reader to the Python
tutorial on the Pickle and cPickle modules.