x = itertools.permutation(range(7))
try:
while True:
print x.next()
except:
pass
x = itertools.permutation(range(7))
for i, e in enumerate(list(x)):
print i+1, e
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)]
x = cycle((1,2,3,4)) while True: print x.next()
tempR = X[:]
for j in range(1, n+1):
minrank = min(tempR)
for k in range(n):
if X[k] == minrank and k< n+1:
X[k] = j
tempR[k] = n+1
return X
4. Name your function, normalsequence(X). .
Here is what I believe a more elegant solution.
def normalsequence(X):
#ensure that the ranks are in sequence!
for i, x in enumerate(sorted(list(set(R[:])))):
for j, y in enumerate(R):
if y == x:
R[j] = i+1
Far more sophisticated (one-liners) for a similar problem is
stack overflow: efficient method to calculate the rank vector of a Python list.
This solution is also part of the Python code for transforming raw scores to ranks in my other post Python, statistics: Ranking a single array of raw scores
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 )
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 )
# -*- 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!
"""
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.
"""
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.070sProblem: Can you rewrite the code so it will be faster? HINT: a comparison is always used inside the loop!
"""
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")
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
myindex = range(11)
for i in range(11):
myindex[i] = order2[order1.index[i]]
# -*- 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
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]
"""
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)
"""
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)) $
def withcommas(s):
p = len(s) -3
s = list(s)
while p > 0:
s.insert(p, ",")
p -= 3
return "".join(s)
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)
(dp1 S'a' S'The quick brown fox' p2 sS'b' S'Hello world.' p3
"""
"""
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")
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.
void showMessage(char *str)
{
puts(str);
}
main()
{
showMessage("Hello World");
}
#includevoid showMessage(char *str) { puts(str); } int main() # explicit declaration. { showMessage("Hello World"); return 0; }
.file "hello.c"
.text
.globl showMessage
.type showMessage, @function
showMessage:
.LFB0:
.cfi_startproc
pushq %rbp
.cfi_def_cfa_offset 16
movq %rsp, %rbp
.cfi_offset 6, -16
.cfi_def_cfa_register 6
subq $16, %rsp
movq %rdi, -8(%rbp)
movq -8(%rbp), %rax
movq %rax, %rdi
call puts
leave
ret
.cfi_endproc
.LFE0:
.size showMessage, .-showMessage
.section .rodata
.LC0:
.string "Hello World"
.text
.globl main
.type main, @function
main:
.LFB1:
.cfi_startproc
pushq %rbp
.cfi_def_cfa_offset 16
movq %rsp, %rbp
.cfi_offset 6, -16
.cfi_def_cfa_register 6
movl $.LC0, %edi
call showMessage
movl $0, %eax
leave
ret
.cfi_endproc
.LFE1:
.size main, .-main
.ident "GCC: (Ubuntu/Linaro 4.4.4-14ubuntu5) 4.4.5"
.section .note.GNU-stack,"",@progbits
The assembler file is one reason most programmers (run away) do not like assembly language! ha ha. Unless of course if the paying job calls for using assmbly (or assembler) langauge. ! But this feature of the gcc to ouput assembler will come in handy when we program a Forth language translator.
<html>
<head>
<title> Hello World Example</title>
<script language="javascript">
function showMessage(message)
{
document.getElementById('location').innerHTML= message;
}
</script>
</head>
<body onload="showMessage('Hello World!')">
<h1>A first Javascript</h1>
A message to one and all:
<span id="location"></span>
</body>
</html>
"""
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.
global _start section .data hello db "Hello, World!", 10 length equ $-hello section .text _start: mov eax, 4 ; write to file mov ebx, 1 ; STDOUT handle mov ecx, hello ; our message mov edx, length ; size of message int 80h ; execute the syscall xor ebx, ebx ; send 0 as 'exit code' mov eax, 1 ; terminate process int 80h ; execute the syscall
section .data hello db "Hello, World!", 10 length equ $-hello
mov eax, 4 ; write to file mov ebx, 1 ; STDOUT handle mov ecx, hello ; our message mov edx, length ; size of message int 80h ; execute the syscall
xor ebx, ebx ; send 0 as 'exit code' mov eax, 1 ; terminate process int 80h ; execute the syscall