473,847 Members | 1,884 Online
Bytes | Software Development & Data Engineering Community
+ Post

Home Posts Topics Members FAQ

numpy help

Howdy,

I'm a college student and for one of we are writing programs to
numerically compute the parameters of antenna arrays. I decided to use
Python to code up my programs. Up to now I haven't had a problem,
however we have a problem set where we are creating a large matrix and
finding it's inverse to solve the problem. To invert the matrix I've
tried using numpy.numarray. linear_algebra. inverse and
numpy.oldnumeri c.linear_algebr a.inverse which both give me the same
error ( I was hoping they called different routines but I think they
call the same one ).
This is the error message I receive:

Traceback (most recent call last):
File "C:\Documen ts and Settings\Chris &
Esther\Desktop\ 636_hw5_2\elen6 36_hw5_2.py", line 60, in <module>
matrix_inverse =
numpy.numarray. linear_algebra. generalized_inv erse(matrix)
File
"C:\Python25\li b\site-packages\numpy\ oldnumeric\line ar_algebra.py", line
59, in generalized_inv erse
return linalg.pinv(a, rcond)
File "C:\Python25\li b\site-packages\numpy\ linalg\linalg.p y", line
557, in pinv
u, s, vt = svd(a, 0)
File "C:\Python25\li b\site-packages\numpy\ linalg\linalg.p y", line
485, in svd
a = _fastCopyAndTra nspose(t, a)
File "C:\Python25\li b\site-packages\numpy\ linalg\linalg.p y", line
107, in _fastCopyAndTra nspose
cast_arrays = cast_arrays + (_fastCT(a.asty pe(type)),)
TypeError: can't convert complex to float; use abs(z)

I've tried inverting small complex matrices and it worked fine. Does
anyone know why it won't work for this larger matrix? Any ideas how I
can work around this problem and get the correct inverse matrix?

Chris

P.S. elen636_math.py is my personal library of functions I've create to
solve the problem while elen636_hw5_2.p y is the program that I'm
actually running
# Purpose:
# This is a library of functions for ELEN 636 that
# so far has the ability to calculate the Sine and
# Cosine integrals as well as the mutual impedance
# between two parallel antennas.
#
# Author: Christopher Smith
# E-mail: mc*****@tamu.ed u
# Date: 10/30/2006

###############
### NOTE: The functions below for the sine and cosine integrals are similar
### to the functions I turned in for homework assignment 4 problem 6
### except that I added the ability to check for convergence.
### I also added the factor into the mutual impedance formula so that the
### answer is given in terms of the terminal input impedance instead of
### the loop impedance as it was formerly giving.
###############

# depends on the math library
from math import *
import numpy.numarray, numpy

def factorial(n):
"""
This function calculates the factorial of a number.
"""
sum = 1.0
for m in range(1, int(n)+1):
sum = float(m)*sum
return sum

def Si(x):
"""
This function computes the sine integral. It uses a power series
expansion that can be found in Abramowitz and Stegun's math
functions reference book.
"""
start = 0.0
stop = 10.0
sine_int = 0.0
convergence = 1.0*10**(-6) # want to have the difference between
# the last run and this run below
# this value
while 1:
for n in range(int(start ), int(stop)):
n = float(n)
sine_int += ((-1)**n)*x**(2*n +1)/((2*n+1)*factor ial(2*n+1))
sine_int_new = sine_int + ((-1.)**stop)*x**( 2.*stop +1.)/((2.*stop+1.)*f actorial(2.*sto p+1.))
converge_check = sine_int_new - sine_int
if abs(converge_ch eck) < convergence:
break
else:
start = stop
stop += 5.0
return sine_int_new

def Ci(x):
"""
This function computes the cosine integral. It uses a power series
expansion that can be found in Abramowitz and Stegun's math
functions reference book.
"""
start = 1.0
stop = 10.0
convergence = 1.0*10.**(-6) # want to have the difference between
# the last run and this run below
# this value
# The first number in the sum is Euler's constant to 10 digits
cosine_int = 0.5772156649 + log(x)
while 1:
for n in range(int(start ), int(stop)):
m = float(n)
cosine_int = cosine_int +((-1)**m)*x**(2*m)/((2*m)*factoria l(2*m))
cosine_int_new = cosine_int + ((-1)**stop)*x**(2 *stop)/((2*stop)*facto rial(2*stop))
converge_check = cosine_int_new - cosine_int
if abs(converge_ch eck) < convergence:
break
else:
start = stop
stop += 5.0
#print stop
return cosine_int_new
def mutual_impedanc e(length1_tot, length2_tot, stagger, d):
"""
This function computes the mutual impedance between two antennas
for the Parallel in Echelon Configuration. The formulas are taken
from a paper by Howard King, "Mutual Impedance of Unequal Length
Antennas in Echelon"

NOTE: all measurements should be entered in wavelengths
"""

# stagger (this is the vertical separation between antenna centers)
# d (this is the horizontal separation between the antennas)

# length1 and length2 are the half length of the antennas, this is
# to conform to King's formulas
length1 = length1_tot/2.0
length2 = length2_tot/2.0

# vertical separation between center of antenna 1 and bottom of antenna 2
h = stagger - length2

# wave propagation constant
beta = 2.0*pi

# formulas to put into mutual impedance equation
u0 = beta*(sqrt(d**2 +(h -length1)**2) +(h -length1))
v0 = beta*(sqrt(d**2 +(h -length1)**2) -(h -length1))
u0prime = beta*(sqrt(d**2 +(h +length1)**2) -(h +length1))
v0prime = beta*(sqrt(d**2 +(h +length1)**2) +(h +length1))
u1 = beta*(sqrt(d**2 +(h -length1 +length2)**2) +(h -length1 +length2))
v1 = beta*(sqrt(d**2 +(h -length1 +length2)**2) -(h -length1 +length2))
u2 = beta*(sqrt(d**2 +(h +length1 +length2)**2) -(h +length1 +length2))
v2 = beta*(sqrt(d**2 +(h +length1 +length2)**2) +(h +length1 +length2))
u3 = beta*(sqrt(d**2 +(h -length1 +2.0*length2)** 2) +(h -length1 +2.0*length2))
v3 = beta*(sqrt(d**2 +(h -length1 +2.0*length2)** 2) -(h -length1 +2.0*length2))
u4 = beta*(sqrt(d**2 +(h +length1 +2.0*length2)** 2) -(h +length1 +2.0*length2))
v4 = beta*(sqrt(d**2 +(h +length1 +2.0*length2)** 2) +(h +length1 +2.0*length2))
w1 = beta*(sqrt(d**2 +h**2) -h)
y1 = beta*(sqrt(d**2 +h**2) +h)
w2 = beta*(sqrt(d**2 +(h +length2)**2) -(h +length2))
y2 = beta*(sqrt(d**2 +(h +length2)**2) +(h +length2))
w3 = beta*(sqrt(d**2 +(h +2.0*length2)** 2) -(h +2.0*length2))
y3 = beta*(sqrt(d**2 +(h +2.0*length2)** 2) +(h +2.0*length2))
#print u0,v0,u0prime,v 0prime,u1,v1,u2 ,v2,u3,v3,u4,v4 ,w1,y1,w2,y2,w3 ,y3

# real part of the mutual impedance between two antennas
R12 = 15*(cos(beta*(l ength1 -h))*(Ci(u0) +Ci(v0) -Ci(u1) -Ci(v1)) \
+sin(beta*(leng th1 -h))*(-Si(u0) +Si(v0) +Si(u1) -Si(v1)) \
+cos(beta*(leng th1 +h))*(Ci(u0prim e) +Ci(v0prime) -Ci(u2) -Ci(v2)) \
+sin(beta*(leng th1 +h))*(-Si(u0prime) +Si(v0prime) +Si(u2) -Si(v2)) \
+cos(beta*(leng th1 -2.0*length2 -h))*(-Ci(u1) -Ci(v1) +Ci(u3) +Ci(v3)) \
+sin(beta*(leng th1 -2.0*length2 -h))*(Si(u1) -Si(v1) -Si(u3) +Si(v3)) \
+cos(beta*(leng th1 +2.0*length2 +h))*(-Ci(u2) -Ci(v2) +Ci(u4) +Ci(v4)) \
+sin(beta*(leng th1 +2.0*length2 +h))*(Si(u2) -Si(v2) -Si(u4) +Si(v4)) \
+2.0*cos(beta*l ength1)*cos(bet a*h)*(-Ci(w1) -Ci(y1) +Ci(w2) +Ci(y2)) \
+2.0*cos(beta*l ength1)*sin(bet a*h)*(Si(w1) -Si(y1) -Si(w2) +Si(y2)) \
+2.0*cos(beta*l ength1)*cos(bet a*(2.0*length2 +h))*(Ci(w2) +Ci(y2) -Ci(w3) -Ci(y3)) \
+2.0*cos(beta*l ength1)*sin(bet a*(2.0*length2 +h))*(-Si(w2) +Si(y2) +Si(w3) -Si(y3)))

# imaginary part of the mutual impedance between two antennas
X12 = 15*(cos(beta*(l ength1 -h))*(-Si(u0) -Si(v0) +Si(u1) +Si(v1)) \
+sin(beta*(leng th1 -h))*(-Ci(u0) +Ci(v0) +Ci(u1) -Ci(v1)) \
+cos(beta*(leng th1 +h))*(-Si(u0prime) -Si(v0prime) +Si(u2) +Si(v2)) \
+sin(beta*(leng th1 +h))*(-Ci(u0prime) +Ci(v0prime) +Ci(u2) -Ci(v2)) \
+cos(beta*(leng th1 -2.0*length2 -h))*(Si(u1) +Si(v1) -Si(u3) -Si(v3)) \
+sin(beta*(leng th1 -2.0*length2 -h))*(Ci(u1) -Ci(v1) -Ci(u3) +Ci(v3)) \
+cos(beta*(leng th1 +2.0*length2 +h))*(Si(u2) +Si(v2) -Si(u4) -Si(v4)) \
+sin(beta*(leng th1 +2.0*length2 +h))*(Ci(u2) -Ci(v2) -Ci(u4) +Ci(v4)) \
+2.0*cos(beta*l ength1)*cos(bet a*h)*(Si(w1) +Si(y1) -Si(w2) -Si(y2)) \
+2.0*cos(beta*l ength1)*sin(bet a*h)*(Ci(w1) -Ci(y1) -Ci(w2) +Ci(y2)) \
+2.0*cos(beta*l ength1)*cos(bet a*(2.0*length2 +h))*(-Si(w2) -Si(y2) +Si(w3) +Si(y3)) \
+2.0*cos(beta*l ength1)*sin(bet a*(2.0*length2 +h))*(-Ci(w2) +Ci(y2) +Ci(w3) -Ci(y3)))

R12_in = R12/(sin(beta*lengt h1)*sin(beta*le ngth2))
X12_in = X12/(sin(beta*lengt h1)*sin(beta*le ngth2))
impedance = (R12_in, X12_in)
return impedance
def top_row_matrix( length1, length2, stagger, stagger_image, radius, m, n):
"""
This function will find the top row of a mutual impdedance matrix
over a ground plane. From the top row we can find the overall matrix
since it is a block Toeplitz matrix
"""

z = [] # list to store our impedance values

# index to step over for the staggering and separation between antennas
stagger_range = range(0, m)
separation_rang e = range(0, n)

# calculate the mutual impedance values for the real planar array
# the first loop gives us the stagger between rows while the second
# loop gives us the separation between dipoles on the same row
for m in stagger_range:
for n in separation_rang e:
h = stagger*m
d = stagger*n
if d == 0: d = radius
trans = mutual_impedanc e(length1, length2, h, d)
z.append(comple x(trans[0],trans[1]))

# Suppose the real antenna array is a plane in x-y at z = 0 and the
# imaginary antenna array is a plane in x-y at z = some spacing
# Since the mutual impedance would include spacing in the
# x,y and z directions we need to condense it down to two spacings
# so that we can input the spacings into our formulas. For the
# spacing between dipoles we will use the hypotenuse formed by the
# x and z coordinates. The y coordinate will be the staggering between
# dipoles.
for m in stagger_range:
for n in separation_rang e:
separation = sqrt((stagger_i mage)**2+(stagg er*n)**2)
trans = mutual_impedanc e(length1, length2, stagger*m, separation)
z.append(comple x(trans[0],trans[1]))

return z
def Toeplitz(a):
"""
This function takes a list in. The list represents
the top row in a Toeplitz matrix. From this row we
can fill the rest of the matrix and then output it.

Note: This function has a dependency on the numpy
library.
"""

# need the length of the list being input
n = len(a)

# fill a matrix the size we need with zeros
matrix = numpy.zeros((n, n), dtype = type(a))

a_new = []
for i in range(-n+1,n):
a_new.append( a[ abs(i) ] )

for i in range(0,n):
matrix[i,:] = a_new[ (n-1) -i : (2*n -1) -i ]

return matrix
def BlockToeplitz(a , M, N):
"""
This function takes a list of N Toeplitz matrices in
and creates a Block Toeplitz matrix from them.

Note: This function has a dependency on the numpy
library.
This function is also specifically for an array
over a groundplane problem. Which means that
the matrix is twice the length and width of a
normal array matrix. The code would need to be
modified to apply to a normal Block-Toeplitz.
"""
# fill a matrix the size we need with zeros
matrix = numpy.zeros((2* M*N,2*M*N), dtype = type(a))

a_new = []
for i in range(-2*N+1,2*N):
a_new.append( a[ abs(i) ] )

for i in range(0, 2*N):
i_new = i*M
for j in range(0, 2*N):
j_new = j*M
matrix[ i_new:i_new+N, j_new:j_new+N ] = a_new[2*N-1-j+i]

return matrix
def h_array_over_gn dplane(M, N):
"""
This function generates the voltage excitation matrix for an
array over a ground plane. The antennas are parallel to the
ground plane which means that according to image theory the
current along the antenna images are in the opposite direction
from those of the real array.

Note: This function has a dependency on the numpy library.
"""

h = []
for n in range(0, 2*M*N):
if n < M*N:
h.append(1.0)
else:
h.append(-1.0)
h = numpy.numarray. array((h), shape = (2*M*N,1))
return h
# Purpose:
# This program will calculate the impedance matrix
# for a planar antenna array over a ground plane.
# This is for ELEN 636 homework 5 problem 2
#
# Author: Christopher Smith
# E-mail: mc*****@tamu.ed u
# Date: 10/30/06

# library so arrays (matrices) can be used
import numpy.numarray, numpy
# library so we can perform linear algebra operations
# like inverse( )
import numpy.numarray. linear_algebra
# import special functions
from elen636_math import mutual_impedanc e, top_row_matrix, Toeplitz
from elen636_math import BlockToeplitz, h_array_over_gn dplane
# math library
from math import *
# length and radius of dipoles (wvlgths)
# horizontal (same in x and y) stagger of matrix (wvlgths)
# vertical separation between array and it's image (wvlgths)
# or twice the vertical separation between array and ground plane
# and the number of elements in the y and x directions m, n
length1 = 0.5
length2 = length1
radius = 0.0025
stagger = 0.55
stag_image = 0.5
M = 7
N = 7

# function to give us back the top row of the m x n matrix
z = top_row_matrix( length1, length2, stagger, stag_image, radius, M, N)

# We now have a list with 2*M*N elements. Using symmetry we can create
# N NxN arrays which correspond to N Toeplitz matrices for our N x N
# Block Toeplitz matrix of 2*M*N x 2*M*N elements. So from those N
# matrices we can completely fill the matrix.
toeplitz_matric es = [] # a list of N Toeplitz matrices
for i in range(0, 2*M*N, N):
toeplitz_matric es.append( Toeplitz( z[i:i +N] ) )

# Now we fill our Block-Toeplitz matrix
matrix = BlockToeplitz(t oeplitz_matrice s, M, N)

# generate the voltage excitations of all the dipoles
h = h_array_over_gn dplane(M, N)

# now we compute Cn by multiply the inverse mutual impedance matrix
# and multiplying it by the h matrix
matrix_inverse = numpy.numarray. linear_algebra. inverse(matrix)
Cn = numpy.numarray. matrixmultiply( matrix_inverse, h)

input_impedance _matrix = h/Cn

#input_impedanc e_element25 = input_impedance _matrix[24]

# creates and opens files to output results
try:
file = open( 'mutual_impedan ce_results.txt' , "w" )
except IOError, message:
print >sys.stderr, "File could not be opened:" , message
sys.exit(1)

print >file, input_impedance _matrix
file.close()



_______________ _______________ _______________ __
Tutor maillist - Tu***@python.or g
http://mail.python.org/mailman/listinfo/tutor
Nov 3 '06 #1
2 3975
Chris Smith wrote:
Howdy,

I'm a college student and for one of we are writing programs to
numerically compute the parameters of antenna arrays. I decided to use
Python to code up my programs. Up to now I haven't had a problem,
however we have a problem set where we are creating a large matrix and
finding it's inverse to solve the problem. To invert the matrix I've
tried using numpy.numarray. linear_algebra. inverse and
numpy.oldnumeri c.linear_algebr a.inverse which both give me the same
error ( I was hoping they called different routines but I think they
call the same one ).
try scipy.linalg.in v etc.

(you are using outdated old NumPy modules with recent array type which sometimes creates type conflicts.

robert

Nov 3 '06 #2
robert wrote:
Chris Smith wrote:
>>Howdy,

I'm a college student and for one of we are writing programs to
numerically compute the parameters of antenna arrays. I decided to use
Python to code up my programs. Up to now I haven't had a problem,
however we have a problem set where we are creating a large matrix and
finding it's inverse to solve the problem. To invert the matrix I've
tried using numpy.numarray. linear_algebra. inverse and
numpy.oldnume ric.linear_alge bra.inverse which both give me the same
error ( I was hoping they called different routines but I think they
call the same one ).


try scipy.linalg.in v etc.

(you are using outdated old NumPy modules with recent array type which sometimes creates type conflicts.

robert
I tried that one also, and it didn't work.

Chris

Nov 3 '06 #3

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

2
3582
by: Justin Lemkul | last post by:
Hello all, I am hoping someone out there will be able to help me. I am trying to install a program that utilizes NumPy. In installing NumPy, I realized that I was lacking Atlas. I ran into the following problems installing Atlas and NumPy, as I realized that NumPy could be installed using the Mac OSX veclib already built in. If anyone has any ideas on how to fix either of these, I would be most grateful. I am fairly new to Python...
1
1532
by: drife | last post by:
Hello, I use the Python Numeric package extensively, and had been an avid user of the "old" scipy. In my view, both pieces of software are truly first rate, and have greatly improved my productivity in the area of scientific analysis. Thus, I was excited to make the transition to the new scipy core (numpy). Unfortunately, I am experiencing a problem that I cannot sort out. I am running Python 2.4.2 on a Debian box (V3.1), using
13
3938
by: Gary Wessle | last post by:
Hi I am trying to install NumPy in my debian/testing linux 2.6.15-1-686. with no numpy for debian/testing, I am left alone, since the experimental version available by debian will result in a dependency nightmares, so after unpacking the downloaded file "numpy-0.9.6.tar.gz"
15
2537
by: greg.landrum | last post by:
After using numeric for almost ten years, I decided to attempt to switch a large codebase (python and C++) to using numpy. Here's are some comments about how that went. - The code to automatically switch python stuff over just kind of works. But it was a 90% solution, I could do the rest by hand. Of course, the problem is that then the code is still using the old numeric API, so it's not a long term solution. Unfortunately, to switch to...
2
2799
by: robert | last post by:
in Gnuplot (Gnuplot.utils) the input array will be converted to a Numeric float array as shown below. When I insert a numpy array into Gnuplot like that below, numbers 7.44 are cast to 7.0 Why is this and what should I do ? Is this bug in numpy or in Numeric? >>m #numpy array array(, , , ..., ,
1
7527
by: Jianzhong Liu | last post by:
Hello, Guys, I have a question about the linear_least_squares in Numpy. My linear_least_squares cannot give me the results. I use Numpy1.0. The newest version. So I checked online and get your guys some examples. I did like this.
3
2211
by: Duncan Smith | last post by:
Hello, Since moving to numpy I've had a few problems with my existing code. It basically revolves around the numpy scalar types. e.g. ------------------------------------------------ array(, ]) Traceback (most recent call last): File "<pyshell#30>", line 1, in -toplevel-
5
6798
by: Marc Oldenhof | last post by:
Hello all, I'm pretty new to Python, but use it a lot lately. I'm getting a crazy error trying to do operations on a string list after importing numpy. Minimal example: Python 2.5.1 (r251:54863, Apr 18 2007, 08:51:08) on win32 Type "help", "copyright", "credits" or "license" for more information.
3
5969
by: Sean Davis | last post by:
I have a set of numpy arrays which I would like to save to a gzip file. Here is an example without gzip: b=numpy.ones(1000000,dtype=numpy.uint8) a=numpy.zeros(1000000,dtype=numpy.uint8) fd = file('test.dat','wb') a.tofile(fd) b.tofile(fd) fd.close()
0
9879
marktang
by: marktang | last post by:
ONU (Optical Network Unit) is one of the key components for providing high-speed Internet services. Its primary function is to act as an endpoint device located at the user's premises. However, people are often confused as to whether an ONU can Work As a Router. In this blog post, we’ll explore What is ONU, What Is Router, ONU & Router’s main usage, and What is the difference between ONU and Router. Let’s take a closer look ! Part I. Meaning of...
0
10978
Oralloy
by: Oralloy | last post by:
Hello folks, I am unable to find appropriate documentation on the type promotion of bit-fields when using the generalised comparison operator "<=>". The problem is that using the GNU compilers, it seems that the internal comparison operator "<=>" tries to promote arguments from unsigned to signed. This is as boiled down as I can make it. Here is my compilation command: g++-12 -std=c++20 -Wnarrowing bit_field.cpp Here is the code in...
0
10643
jinu1996
by: jinu1996 | last post by:
In today's digital age, having a compelling online presence is paramount for businesses aiming to thrive in a competitive landscape. At the heart of this digital strategy lies an intricately woven tapestry of website design and digital marketing. It's not merely about having a website; it's about crafting an immersive digital experience that captivates audiences and drives business growth. The Art of Business Website Design Your website is...
1
10705
by: Hystou | last post by:
Overview: Windows 11 and 10 have less user interface control over operating system update behaviour than previous versions of Windows. In Windows 11 and 10, there is no way to turn off the Windows Update option using the Control Panel or Settings app; it automatically checks for updates and installs any it finds, whether you like it or not. For most users, this new feature is actually very convenient. If you want to control the update process,...
0
9477
agi2029
by: agi2029 | last post by:
Let's talk about the concept of autonomous AI software engineers and no-code agents. These AIs are designed to manage the entire lifecycle of a software development project—planning, coding, testing, and deployment—without human intervention. Imagine an AI that can take a project description, break it down, write the code, debug it, and then launch it, all on its own.... Now, this would greatly impact the work of software developers. The idea...
1
7879
isladogs
by: isladogs | last post by:
The next Access Europe User Group meeting will be on Wednesday 1 May 2024 starting at 18:00 UK time (6PM UTC+1) and finishing by 19:30 (7.30PM). In this session, we are pleased to welcome a new presenter, Adolph Dupré who will be discussing some powerful techniques for using class modules. He will explain when you may want to use classes instead of User Defined Types (UDT). For example, to manage the data in unbound forms. Adolph will...
0
5907
by: adsilva | last post by:
A Windows Forms form does not have the event Unload, like VB6. What one acts like?
1
4521
by: 6302768590 | last post by:
Hai team i want code for transfer the data from one system to another through IP address by using C# our system has to for every 5mins then we have to update the data what the data is updated we have to send another system
3
3158
bsmnconsultancy
by: bsmnconsultancy | last post by:
In today's digital era, a well-designed website is crucial for businesses looking to succeed. Whether you're a small business owner or a large corporation in Toronto, having a strong online presence can significantly impact your brand's success. BSMN Consultancy, a leader in Website Development in Toronto offers valuable insights into creating effective websites that not only look great but also perform exceptionally well. In this comprehensive...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.