Python, finding nearest matching color name for RGB values

There are 147 standard color names for use in HTML. here is a simple Python code to find the nearest matching color name given input R,G, B values (from 0-255). The minimization criterion is the sum of absolute vaues of the differences from the color r, g, b values of the colornames.

"""
File    colornames.py
Author  Ernesto P. Adorio, Ph.D
Version 0.0.1 February 1, 2012
"""


# Src: http://www.w3schools.com/html/html_colornames.asp
lines = """AliceBlue  #F0F8FF   Shades Mix
AntiqueWhite  #FAEBD7   Shades Mix
Aqua  #00FFFF   Shades Mix
Aquamarine  #7FFFD4   Shades Mix
Azure  #F0FFFF   Shades Mix
Beige  #F5F5DC   Shades Mix
Bisque  #FFE4C4   Shades Mix
Black  #000000   Shades Mix
BlanchedAlmond  #FFEBCD   Shades Mix
Blue  #0000FF   Shades Mix
BlueViolet  #8A2BE2   Shades Mix
Brown  #A52A2A   Shades Mix
BurlyWood  #DEB887   Shades Mix
CadetBlue  #5F9EA0   Shades Mix
Chartreuse  #7FFF00   Shades Mix
Chocolate  #D2691E   Shades Mix
Coral  #FF7F50   Shades Mix
CornflowerBlue  #6495ED   Shades Mix
Cornsilk  #FFF8DC   Shades Mix
Crimson  #DC143C   Shades Mix
Cyan  #00FFFF   Shades Mix
DarkBlue  #00008B   Shades Mix
DarkCyan  #008B8B   Shades Mix
DarkGoldenRod  #B8860B   Shades Mix
DarkGray  #A9A9A9   Shades Mix
DarkGrey  #A9A9A9   Shades Mix
DarkGreen  #006400   Shades Mix
DarkKhaki  #BDB76B   Shades Mix
DarkMagenta  #8B008B   Shades Mix
DarkOliveGreen  #556B2F   Shades Mix
Darkorange  #FF8C00   Shades Mix
DarkOrchid  #9932CC   Shades Mix
DarkRed  #8B0000   Shades Mix
DarkSalmon  #E9967A   Shades Mix
DarkSeaGreen  #8FBC8F   Shades Mix
DarkSlateBlue  #483D8B   Shades Mix
DarkSlateGray  #2F4F4F   Shades Mix
DarkSlateGrey  #2F4F4F   Shades Mix
DarkTurquoise  #00CED1   Shades Mix
DarkViolet  #9400D3   Shades Mix
DeepPink  #FF1493   Shades Mix
DeepSkyBlue  #00BFFF   Shades Mix
DimGray  #696969   Shades Mix
DimGrey  #696969   Shades Mix
DodgerBlue  #1E90FF   Shades Mix
FireBrick  #B22222   Shades Mix
FloralWhite  #FFFAF0   Shades Mix
ForestGreen  #228B22   Shades Mix
Fuchsia  #FF00FF   Shades Mix
Gainsboro  #DCDCDC   Shades Mix
GhostWhite  #F8F8FF   Shades Mix
Gold  #FFD700   Shades Mix
GoldenRod  #DAA520   Shades Mix
Gray  #808080   Shades Mix
Grey  #808080   Shades Mix
Green  #008000   Shades Mix
GreenYellow  #ADFF2F   Shades Mix
HoneyDew  #F0FFF0   Shades Mix
HotPink  #FF69B4   Shades Mix
IndianRed   #CD5C5C   Shades Mix
Indigo   #4B0082   Shades Mix
Ivory  #FFFFF0   Shades Mix
Khaki  #F0E68C   Shades Mix
Lavender  #E6E6FA   Shades Mix
LavenderBlush  #FFF0F5   Shades Mix
LawnGreen  #7CFC00   Shades Mix
LemonChiffon  #FFFACD   Shades Mix
LightBlue  #ADD8E6   Shades Mix
LightCoral  #F08080   Shades Mix
LightCyan  #E0FFFF   Shades Mix
LightGoldenRodYellow  #FAFAD2   Shades Mix
LightGray  #D3D3D3   Shades Mix
LightGrey  #D3D3D3   Shades Mix
LightGreen  #90EE90   Shades Mix
LightPink  #FFB6C1   Shades Mix
LightSalmon  #FFA07A   Shades Mix
LightSeaGreen  #20B2AA   Shades Mix
LightSkyBlue  #87CEFA   Shades Mix
LightSlateGray  #778899   Shades Mix
LightSlateGrey  #778899   Shades Mix
LightSteelBlue  #B0C4DE   Shades Mix
LightYellow  #FFFFE0   Shades Mix
Lime  #00FF00   Shades Mix
LimeGreen  #32CD32   Shades Mix
Linen  #FAF0E6   Shades Mix
Magenta  #FF00FF   Shades Mix
Maroon  #800000   Shades Mix
MediumAquaMarine  #66CDAA   Shades Mix
MediumBlue  #0000CD   Shades Mix
MediumOrchid  #BA55D3   Shades Mix
MediumPurple  #9370D8   Shades Mix
MediumSeaGreen  #3CB371   Shades Mix
MediumSlateBlue  #7B68EE   Shades Mix
MediumSpringGreen  #00FA9A   Shades Mix
MediumTurquoise  #48D1CC   Shades Mix
MediumVioletRed  #C71585   Shades Mix
MidnightBlue  #191970   Shades Mix
MintCream  #F5FFFA   Shades Mix
MistyRose  #FFE4E1   Shades Mix
Moccasin  #FFE4B5   Shades Mix
NavajoWhite  #FFDEAD   Shades Mix
Navy  #000080   Shades Mix
OldLace  #FDF5E6   Shades Mix
Olive  #808000   Shades Mix
OliveDrab  #6B8E23   Shades Mix
Orange  #FFA500   Shades Mix
OrangeRed  #FF4500   Shades Mix
Orchid  #DA70D6   Shades Mix
PaleGoldenRod  #EEE8AA   Shades Mix
PaleGreen  #98FB98   Shades Mix
PaleTurquoise  #AFEEEE   Shades Mix
PaleVioletRed  #D87093   Shades Mix
PapayaWhip  #FFEFD5   Shades Mix
PeachPuff  #FFDAB9   Shades Mix
Peru  #CD853F   Shades Mix
Pink  #FFC0CB   Shades Mix
Plum  #DDA0DD   Shades Mix
PowderBlue  #B0E0E6   Shades Mix
Purple  #800080   Shades Mix
Red  #FF0000   Shades Mix
RosyBrown  #BC8F8F   Shades Mix
RoyalBlue  #4169E1   Shades Mix
SaddleBrown  #8B4513   Shades Mix
Salmon  #FA8072   Shades Mix
SandyBrown  #F4A460   Shades Mix
SeaGreen  #2E8B57   Shades Mix
SeaShell  #FFF5EE   Shades Mix
Sienna  #A0522D   Shades Mix
Silver  #C0C0C0   Shades Mix
SkyBlue  #87CEEB   Shades Mix
SlateBlue  #6A5ACD   Shades Mix
SlateGray  #708090   Shades Mix
SlateGrey  #708090   Shades Mix
Snow  #FFFAFA   Shades Mix
SpringGreen  #00FF7F   Shades Mix
SteelBlue  #4682B4   Shades Mix
Tan  #D2B48C   Shades Mix
Teal  #008080   Shades Mix
Thistle  #D8BFD8   Shades Mix
Tomato  #FF6347   Shades Mix
Turquoise  #40E0D0   Shades Mix
Violet  #EE82EE   Shades Mix
Wheat  #F5DEB3   Shades Mix
White  #FFFFFF   Shades Mix
WhiteSmoke  #F5F5F5   Shades Mix
Yellow  #FFFF00   Shades Mix
YellowGreen  #9ACD32"""

def rgbfromstr(s):
    # s starts with a #.
    r, g, b = int(s[1:3],16), int(s[3:5], 16),int(s[5:7], 16)
    return r, g, b

def findnearestcolorname(R,G, B, colorD):
    mindiff = None

    for d in colorD:
       r, g, b = rgbfromstr(colorD[d])
       diff = abs(R -r)*256 + abs(G-g)* 256 + abs(B- b)* 256 
       if mindiff is None or diff < mindiff:
          mindiff = diff
          mincolorname = d
    return mincolorname
 

if __name__ == "__main__":
  D={}
  for line in lines.split("\n"):
    tokens = line.split()
    if len(tokens) > 1:
        D[tokens[0]]= tokens[1]
     
  R, G, B = 0, 127, 127
  colorname = findnearestcolorname(0, 127, 127, D)
  print hex(R), hex(G), hex(B)
  print colorname, D[colorname]
 

When the above program runs, it outputs


$ python colornames.py
0x0 0x7f 0x7f
Teal #008080



Questions.

1. Modify the minimizaiton criterion to use least sum of squares of the differences.

2. Provide error checking for input R, G, B values to be within the interval [0, 255].

3. Procide different ways to specify the input RGB values, aside from the (R, G, B) values as done in the program above.

Python: Permutations, Combinations, Cartesian Products with Itertools module

Python itertools module has built-in functions for permutations, combinations, cycles and products.
So the Python programmers can focus on building the next killer-apps :) which will dominate the whole world, instead
of reading obtuse treatises, papers, and other publications on algorithms for generating combinatorial objects.

Here we list the combinatorial capabilities of Python (versions 2.6+) in the itertools module.

• Permutations

  
  itertools.permutation create a permutation object from an iterable input or string. You can list
  the successive permutations by using object.next() or by creating a list (a bad idea if the size of the object is large!)
  list(object). Here is an example:
  
  

  x = itertools.permutation(range(7))
     try:
        while True:
          print x.next()
     except:
        pass
  

  
  You can also do:
  
  

  x = itertools.permutation(range(7))
     for i, e in enumerate(list(x)):
         print i+1, e
  

  
  Converting to a list give rise to a dangerous situation where you can run out of memory!. but if you know that it can be safely stored,
  by all means use a list conversion.
  
  Permutations can also accept an extra parameter r so you can generate all permutations of size n of the object taken r at a time.
  The number of elements in such permutations is given by n!/(n-r)!. If r is not specified, then the number of permutations is n!
  which look innocous.(50! is nothing to sneeze at).
  
  

• combinations

itertools.combinations creates a combinatorial object. Here ordering does not matter. As an example, the combinations of 6 elements taken 3 at a time can be listed in two groups, as follows:

import itertools

seq = range(6)
x = itertools.combinations(seq, 3)

try:
  for i, e in list(enumerate(x)):
     f = [y for y in seq if y not in e]
     print i+1, list(e), f
except:
   pass

The output of permutations and combinations are tuples for fast execution of the functions, and have to be typecast to a list. Itertools also provide combinations with replacement. The best way to see the difference with the standard combinations is to see them in action from the command line:

x = combinations_with_replacement([1,2,3,4], 2)
y = combinations([1,2,3,4], 2)

list(x)
[(1, 1),
 (1, 2),
 (1, 3),
 (1, 4),
 (2, 2),
 (2, 3),
 (2, 4),
 (3, 3),
 (3, 4),
 (4, 4)]

list(y)
[(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]

• Cartesion products

  
  Itertools has cartesian products! So instead of writing (this from the documentation)
  
  ((x,y) for x in A for y in B)
  
  we simply write product(A,B)
  
  Wait, how about the cartesian product of A with itself?
  you can write product(A,A) or if you want more, product(A,A,A,A). The latter can be written as
  product(A, repeat=4). Python programmers can be lazy as Perl programmers!
  
  

• cycle





Not strictly referrring to a mathematical object. Rather it refers to the way an element in an
iterable is accessed. Instead of stopping at the last element, the access is cycled back to the first element.
Try this on the command line (not running in a script file)



x = cycle((1,2,3,4))
while True:
   print x.next()


You will notice that this loop never terminates, and after the 4 is printed, the next element to be printed is the first element.
Type Ctrl-C to abort the printing.

Aside from these basic combinatorial objects, Itertools provide now a repeat function! For example, a list of 100
elements with all values equal to 5: [5] * 100 can now be written as itertools.repeat(5, 100). Easier to understand for beginners.


Dont you just love Python?

Hopefully, Python will offer genuine cycle objects (where you can multiply two cycles) and derangements with basic functions such as nPr and nCr. These combinatorial features are solutions waiting for problems and applications!

TO BE CONTINUED.

Python: Computing sunrise sunset times with ephem

It is incredibly easy to find the sunrise and sunset for the city of Manila using Python and ephem libary. Here is for today, September 28, 2011

import ephem

manila = ephem.city("Manila")
sun    = ephem.Sun()
manila.date="2011/09/28"
sunrise= manila.next_rising(sun)
sunset = manila.next_setting(sun)

print ephem.localtime(sunrise)
print ephem.localtime(sunset )

When run within ipython, it results in


In[xxx]: print ephem.localtime(sunrise)
2011-09-29 05:45:16.000003

In [xxx]: print ephem.localtime(sunset )
2011-09-28 17:48:28.000002



Times returned by ephem are always in UTC, a point which puzzled me for more than 30 minutes, wondering what was wrong! We only have to explicitly call the localtime function.

The ephem library knows some cities already. Here is information about Manila:


ephem.Observer date='2011/9/28 00:00:00' epoch='2000/1/1 12:00:00' long=120:58:56.0 lat=14:36:15.0

elevation=7.9248m horizon=0:00:00.0 temp=15.0C pressure=1012.29834492mBar



What if you live in Angeles City, the Philippines? Here is a modified Python code:

import ephem

# manila = ephem.city("Manila")
# Angeles has longitude, latitude of 15°09'N	120°33'E respectively.
mycity= ephem.Observer()
mycity.long= "120:33"
mycity.lat="15:09"

sun    = ephem.Sun()
mycity.date="2011/09/28"
sunrise= mycity.next_rising(sun)
sunset =mycity.next_setting(sun)
print ephem.localtime(sunrise)
print ephem.localtime(sunset )

Here is the output for Angeles City.


In [xxx]: print ephem.localtime(sunrise)
2011-09-29 05:47:04.000003

In [xxx]: print ephem.localtime(sunset )
2011-09-28 17:50:07.000002


We shall put up an ephemeris data page shortly. To view instructions for installation of PyEphem, please visit
http://free-software-explorations.blogspot.com/2010/07/pyephem-ephemeris-computation-library.html

Python: basic conjugate gradient linear solver

Visit Wikipedia:Conjugate gradient method for a clear and succinct explanation of the conjugate gradient method for solving linear equations. We based our Python code on the algorithm described in this article. We are pleased that we were able to get the same results for the simple example involving a 2 by 2 system,

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

"""
File     conjgrad.py
Author   Ernesto P. Adorio
         UPDEPP Clarfield
email    ernesto.adorio@gmail.com          
ref.     http://en.wikipedia.org/wiki/Conjugate_gradient_method
"""

from matlib import *

def conjgrad(A, b, x = None, ztol = 1.0e-4, maxiter = None):
    """
    Solves the equation Ax = b by iteration using conjugate gradients.
    Arguments 
      A- symmetric positive definite matrix.

    """
    print "input A=", A
    print "input b=", b
    print "input x=", x
    if x is None:
       x = [0] * len(b)
    m, n = matdim(A)
    if m != n or m!= len(b):
       raise ValueError
    if maxiter is None:
       maxiter = m    
    rk = vecsub(b, matvec(A, x))
    print "r_0 =", rk
    pk = rk[:]
    num    = dot(rk, rk)

    for k in range(1,maxiter):
        print "iteration #", k
        Apk    = matvec(A, pk)
        den    = dot(pk, Apk)
        alphak = num/float(den)
        print "num/den=", num, den, alphak

        #update x.
        x  = [e  + alphak * p for (e,p) in zip(x, pk)]
        print "x = ", x

        #update rk.
        rk    = [r - alphak * apk for r, apk in zip(rk, Apk)]
        dotrk = dot(rk, rk)
        if num < ztol:
           return x, k
        print "rkp1=", rk

        #update pk
        betak = dotrk/num
        num   = dotrk
        pk    = [r + betak* p for r, p in zip(rk, pk)] 
        print "betak=", betak
        print "pk=", pk
        print
    return x, maxiter  

def vec2mat(x, nrows, ncols):
    if nrows * ncols != len(x):
       raise ValueError,"incompatible dimensions in vec2mat()"
    M = []
    for i in range(nrows):
       M.append(x[i*ncols:(i+1)*ncols])

    return M   

     
if __name__ == "__main__":
   A = vec2mat([4,1,1,3],2,2)
   b = [1,2]
   x = [2,1]
   sol = conjgrad(A,b, x,ztol = 1.0e-4, maxiter=5)
   print "sol=", sol

We will soon add the vec2mat() function included in the code above in the matlib.py function. Please note that this author recommends using the numpy and scipy Python modules. They are well tested and provide a wrapper to the linear algebra routines developed in EISPACK and LINPACK. scientific libraries which were coded in C. Here are the results of running the above from the terminal:

toto@toto-Aspire-4520:~/Blogs/lm$ python conjgrad.py 
input A= [[4, 1], [1, 3]]
input b= [1, 2]
input x= [2, 1]
r_0 = [-8, -3]
iteration # 1
num/den= 73 331 0.220543806647
x =  [0.2356495468277946, 0.33836858006042303]
rkp1= [-0.2809667673716012, 0.7492447129909365]
betak= 0.00877136937414
pk= [-0.3511377223647101, 0.7229306048685207]

iteration # 2
num/den= 0.640309964312 1.55338036659 0.412204234122
x =  [0.09090909090909094, 0.6363636363636365]
rkp1= [5.551115123125783e-17, 0.0]
betak= 4.81249407751e-33
pk= [5.551115123125783e-17, 3.4790992543769425e-33]

iteration # 3
num/den= 3.08148791102e-33 1.23259516441e-32 0.25
x =  [0.09090909090909095, 0.6363636363636365]
sol= ([0.09090909090909095, 0.6363636363636365], 3)
toto@toto-Aspire-4520:~/Blogs/lm$         
 

Of course remove the print statements to avoid the computational progress reports. Our aim here is mainly pedagogical. The Python code is easier to read and one can use a debugger (as in Eric4) to walk through the code!

Does a string represents a floating point number, an integer or a plain string?

I have been programming in Python a long long time, I have not coded as a function one of the most common tests on a string, for the simple reason that it is easy to program directly. Here, we ask: does a string represents a floating point number(reals), an integer or if not these types, a plain string?
Here are simple test cases: -123., -123e0, -123, . will be considered reals or floats. Numbers which do not contain a decimal point or dot or an exponentiation part will be considered integers. Here is our function.

"""
file            dtypeval.py

author      Ernesto P. Adorio, PhD.
                 UPDEPP at Clarkfield, Pampanga

version    0.0.1 september 19, 2011
"""



def dtypeval(s, numericfloat=False):
    """
    Simple function to return the data type and value of a string.
    """
    try:
        if numericfloat or ("." in s) or ("e" in s) or ("E" in s):
           f = float(s)
           return "float", f
        return "int", int(s) 
    except:
        # print "[%s] is not a number!" % s
        return "str", s


if __name__ == "__main__":
   teststring  = ["-123.", "-123", "1e5", "1E5", ".001", "ABC"]

   for s in teststring:
       print s,"==>", dtypeval(s)

Save the above to dtypeval.py and run from the command line python dtype.val. Here is the output:

-123. ==> ('float', -123.0)
-123 ==> ('int', -123)
1e5 ==> ('float', 100000.0)
1E5 ==> ('float', 100000.0)
.001 ==> ('float', 0.001)
ABC ==> [ABC] is not a number!
('str', 'ABC')

Note that numericfloat is set to true if the function is to return float always even if the input sting represents an integer.

Python: Sampling from a normal population

Here is code to create N samples of size n from a normal population with mean mu and standard deviation sigma.

"""
File             normalsim.py
Author           Dr. Ernesto P. Adorio
Desc             Sampling from the normal distribution.
Version          0.0.1 Sep. 1, 2011
"""



import numpy as np
import scipy.stats as stat
from math import sqrt


def  normalsim(mu, sigma, n, N, statistic="mean"):
    out = [0] * N
    for i in range(N):
        sample = stat.norm.rvs(mu, sigma, size=n)
        #print sample, 
        if statistic=="mean":
           out[i] = np.mean(sample)
        elif statistic=="svar":
           out[i] = np.var(sample, ddof=1)
        elif statistic=="pvar":
           out[i] = np.var(sample)
        elif statistic=="ssdev":
           out[i] = sqrt(np.var(sample, ddof=1))
        elif statistic=="psdev":
           out[i] = np.var(sample)
        elif statistic== "min":
           out[i] = min(sample)
        elif statistic == "max":
           out[i] = max(sample)
        else:
           raise ArgumentError,"bad statistic."
        #print out[i]
    return out

v = normalsim(10, 4, 10, 1000, statistic="svar")
print np.mean(v), np.var(v)

When I run the program under the time command, time python normalsim.py, the shell outputted

15.7770857263 55.2458134956

real    0m2.681s
user    0m1.080s
sys     0m0.070s

Problem: Can you rewrite the code so it will be faster? HINT: a comparison is always used inside the loop!

Python: converting names to surname, firstname format.

we want to convert names like

Jose P. Rizal to Rizal, Jose P.
Jose P. Rizal Jr. to Rizal Jr.,   Jose P.
Ma. Jose P. Rizal III to Rizal III, Ma. Jose P.

Here is Python code to perform the conversion. I am wondering if one can write a more elegant version, where by elegant, it does the conversion in a faster, shorter way!

"""
file     lastfirstname.py
author   Ernesto P. Adorio, PhD.
desc     Converts name to last name, first name format.
version  0.0.1 july 27, 2011
"""

def makestr(sarray, lastindex):
    """
   
    """
    output = ""
    for s in sarray[:lastindex]:
       output += (" " + s )
    return output

def surnamefirstname(name):
    """
    Translates a name in firstname surname to
          surname, firstname format.
    Complications arise because surname can include "Jr." or romanized numbers "III"
    """
    parts = name.split()
    n = len(parts)
    if n == 2:
       # easiest case
       return parts[1] + ", " + parts[0]
    elif n > 2:
       if parts[-1].lower() in ["jr.", "jr"]:
          return parts[-2] +" Jr.," + makestr(parts, -2)
       elif parts[-1].lower() in ["ii", "iii", "iv", "v"]:
          return parts[-2] + " " +parts[-1].upper()+ ","  + makestr(parts, -2)
       else:
          return parts[-1] + "," + makestr(parts, -1)
 
print surnamefirstname("Ma. Jose P. Rizal III")

 
       

Post you own better version in the comments. We shall learn together to improve our Python skills.

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.

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.

Python: a high level word count program

Introduction

Python allows the beginning programmer to use powerful built-in powerful routines, freeing the programmer from imlementing on his own (usually error prone versions)! Here, we use Python A dictionarie, convenient string splitting methods to do the non-trivial word counting problem. Non-trivial in the sense it does more than printing a message "Hello World" to the screen.

the Python code

"""
"""
file     wordcount.py
author   Ernesto P. Adorio
         U.P.Clarkfield
version  0.0.1
desc     counts words, lines
         lines are setoff by new line character, "\n"
         quoted strings are stripped off of quotes.   
"""

def wc(textfilename):
    D = {}
    nc = 0 # number of characters.
    nw = 0 # number of words.
    nl = 0 # number of lines.
    for line in open(textfilename, "r").readlines():
        # print line #decomment for debugging.
        nc += len(line)
        words = line.split()
        for word in words:
            if word in D:
               D[word] += 1
            else:
               D[word] = 1
            nw += 1
        nl += 1
    return nc, nw, nl


if __name__ == "__main__":
   print wc("wordcount.py")

Notice the code's compactness. In a lower level language such as C, the code will be far longer and complicated. The essence of the code is:

read in all lines and store in an array.

process each line.
add length of line to number of characters.
split the line into words, which are tokens separated by whitespaces.
add number of words in line the number of total words.
update the line number.

When the program is fed the code for itself, it printed the number of c
haracters, words and lines as


$ python wordcount.py
(762, 104, 36)


As a check the wc utility produces ( wc wordcount.py) in reverse order


36 104 762

Issues?

There may be instances when the number of lines is off by one. This happens if the last line does not contains a terminating line feed character.
Another thing is that the code may be feed an extremely large file. If the working memory is small, it may run out of space! This problem is easy to fix and we will revise the above initial working code later.

Python: Converting a matrix of raw scores to ranks.

There are instances where we have to convert a matrix of raw scores to ranking scores, for example in Statistics where we have to compute the Kruskal-Wallis statistics of a block of values. Although we have written a Kruskal Wallis routine previously, we shall write an independent function to return ranks based on raw scores.

A complication in computing ranks is that of processing ties. An approach is to replace any tied scores with the mean of all ranks. For example if 4 values v1,v2,v3,v4 have ranks 5,6,7,8, then each value should have an associated rank of (5+6+7+8)/4.

Here is an outline for converting raw scores to ranks:
1. Create a list L containing tuples of the form (v, r,c) for value, row, column for each element in the input raw scores matrix.
2. Sort the temporaty list L in ascending order.
3. Process the list L to fill in the output Ranks matrix,

Here is our python code for outputting ranks values, complete with an example from Berensonś Basic Business Statistics, page 493.(We do have a local copy of the book) for testing the routine.

"""
File    matrix_scores_to_ranks.py
Author  Dr. Ernesto P. Adorio
        U.P. Clarkfield
        ernesto.adorio@gmail.com
"""


def matrix_scores_to_ranks(M, break_ties = 1, start_rank = 1, ZTOL = 1.0e-5):
    """
    Args
       M          - a list of lists of values.
       break_ties - 1 if ties need special processing, otherwise zero
       start_rank - 0 or 1
       ZTOL       - equality comparison of two real numbers.
    Return value
       A list of lists of ranks of the values. 
    """
    R=[]
    L=[]
    for i, row in enumerate(M):
        R.append([0] * len(row))
        for j, v in enumerate(row):
    L.append((v, i, j))
    L.sort()

    if break_ties != 1:
       for k, tup in enumerate(L):
          v, i, j = tup
   if break_ties != 1:   
     R[i][j] = k + start_rank
    else:
        prevv = L[0][0]    
        starti = 0 
        endi   = 0 
        for k, tup in enumerate(L):
            v, i, j = tup
            if abs(v- prevv) < ZTOL: # a tie?
               endi = k
            else:
               # print previous tied scores!
               rank = sum(range(starti, endi+1))/(endi-starti + 1.0) + start_rank

               for i in range(starti, endi+1):
                   ii = L[i][1]
                   jj = L[i][2]
                   R[ii][jj] = rank
               prevv = v
               starti = endi = k   
        # Final processing of any remaining ties.
        rank = sum(range(starti, endi+1))/(endi-starti + 1.0) + start_rank
        for i in range(starti, endi+1):
            ii = L[i][1]
            jj = L[i][2]
            R[ii][jj] = rank
    return R

#Berenson p. 493.

if __name__ == "__main__":
  M = [[ 18.5, 24.0, 17.2, 19.9, 18.0],
     [ 26.3, 25.3, 24.0, 21.2, 24.5],
     [20.6,  25.2, 20.8, 24.7, 22.9],
     [25.4,  19.9, 22.6, 17.5, 20.4]]

  R = scores_to_ranks(M, break_ties = 1, start_rank = 1)

  for row in R:
     print row

When the above program runs, it outputs and matches with the results in the textbook.

[4.0, 13.5, 1.0, 5.5, 3.0] [20.0, 18.0, 13.5, 10.0, 15.0] [8.0, 17.0, 9.0, 16.0, 12.0] [19.0, 5.5, 11.0, 2.0, 7.0]

The above is for a matrix. We will post array_scores_to_ranks() with the next article posting. Of course, the output must be transposed to read properly. start_rank may be set to any starting value, usually a 1, but there may be occasion when a starting zero rank may be appropriate. The above code is written for clarity. We may speed it up later, but there is no incentive to do so. We think most of the running time is spent on the sort, and Python sort routine is already fast enough! Feb. 18, 2011: I need to install the Syntax Highlighter addin immediately! Done. I refer a blogger user to Coding Freak There are other strategies in breaking ties. See Python, statistics: Ranking a single array of raw scores and we will update the above code someday.