I'm trying to create a piece of software that calculates the density distribution inside a white dwarf by dividing it into a bunch of layers and running a crude physics simulation on them. I wrote some functions that I tried to vectorize, but I keep getting the error TypeError: return arrays must be of ArrayType.
From what I've been able to gather from other questions on this site, the culprit is probably running operations on entities with different types. So I tried rewriting the frompyfunc() statements to vectorize statements and tried setting the type to be np.double, "double", and "float64", but I ended up getting "invalid otype." And yes, I remembered to get rid of the ninput and noutput arguments when doing this.
I also tried setting all of my arrays to the type double and clarifying which variables were global in my original functions, but neither of those helped either.
This is the function
def selfBelow(m,r,x): #The gravitational self-force of a layer below the border
if x>r:
return 3*G*m**2*x*(5*r**3 + 6*r**2*x + 3*r*x**2 + x*3) / (5*(r**2 + r*x + x**2)**3)
else:
return 1/0 #Return an error in this case because we need x to be greater than r
This is my attempt to vectorize it: FSB = np.frompyfunc(selfBelow,2,1)
And this is my attempt to call the vectorized function: forces[0]-=FSB(masses[0], 0, radii[0])
When I run it, I get this error:
Traceback (most recent call last):
File "C:\...\White Dwarf.py", line 78, in <module>
forces[0]-=FSB(masses[0], 0, radii[0])
TypeError: return arrays must be of ArrayType
Here's the full code in case you need it:
import numpy as np
import matplotlib as plt
h=1.054571817e-34 #Reduced Planck constant
G=6.67430e-11 #Gravitational constant
m_e=9.10938370153e-31 #Mass of an electron in kg
u=1.66053907e-27 #Atomic mass unit in kg
m_Sun=1.9885e30
r_Earth=6.3781e6
n=2 #Number of baryons per electron
C=100 #Number of layers
def volume(r, R):
return 4*np.pi/3*(R**3-r**3)
def selfBelow(m,r,x): #The gravitational self-force of a layer below the border
if x>r:
return 3*G*m**2*x*(5*r**3 + 6*r**2*x + 3*r*x**2 + x*3) / (5*(r**2 + r*x + x**2)**3)
else:
return 1/0
def selfAbove(m,x,R): #The gravitational self-force of a layer abpve the border
return 9*G*m**2*x**2*(R**2 + 3*R*x + x**2)/(10*(R**2+R*x+x**2)**3)
def elseBelow(M,m,r,x): #The gravitational force energy from the cumulative masses on the layer below
if x>r:
return 3*G*M*m*x*(2*r+x)/(2*(r**2+r*x+x**2)**2)
else:
return 1/0
def elseAbove(M,m,x,R): #The gravitational force from the cumulative masses on the layer below
if x>r:
return 3*G*M*m*x*(2*R+x)/(2*(R**2+R*x+x**2)**2)
else:
return 1/0
def fermiBelow(n,r,x): #The fermi pressure from the layer below
if x>r:
return -(9*pow(3/2,1/3)*pow(np.pi,2/3)*h**2*x**2)/(10*m_e)*pow(n/(x**3-r**3),5/3)
else:
return 1/0
def fermiAbove(n,x,R):
if R>x:
return -(9*pow(3/2,1/3)*pow(np.pi,2/3)*h**2*x**2)/(10*m_e)*pow(n/(x**3-r**3),5/3)
else:
return 1/0
V = np.frompyfunc(volume,2,1)
FSB = np.frompyfunc(selfBelow,2,1)
FSA = np.frompyfunc(selfAbove,2,1)
FEB = np.frompyfunc(elseBelow,2,1)
FEA = np.frompyfunc(elseAbove,2,1)
FFB = np.frompyfunc(fermiBelow,2,1)
FFA = np.frompyfunc(fermiAbove,2,1)
sunMasses=1
M=m_Sun*sunMasses
Particles=M/(n*u)
startRadius = h**2/(G*M**2*m_e)*pow(Particles,5/3)*pow(9*np.pi/4,2/3)
radii = startRadius*np.arange(1,C+1)/C
masses = np.ndarray(C, np.double)
masses[0:C] = M*V(np.arange(0,C),np.arange(1,C+1))/volume(0,C)
cumuMasses = np.ndarray(C, np.double)
cumuMasses[0:C] = M*V(0,np.arange(1,C+1))/volume(0,C) #Cumulative mass of all layers below the border
N = np.ndarray(C,np.double)
N[0:C] = Particles*V(np.arange(0,C),np.arange(1,C+1))/volume(0,C)
forces =np.zeros(C, np.double)
print(radii.dtype, masses.dtype, cumuMasses.dtype, N.dtype, forces.dtype)
for i in range(C):
print(i)
forces.fill(0)
forces[0]-=FSB(masses[0], 0, radii[0]) #Special case because radii does not include 0
forces[1:C]-=FSB(masses[1:C], radii[0:-1], radii[1:C]) #The 0:-1 index means up to, but not including, the last element
#0 is not included for reason above
forces[0:-1]-=FSA(masses[1:C], radii[0:-1], radii[1:C])
forces[1:C]-=FEB(cumuMasses[0:-1], masses[1:C], radii[0:-1], radii[1,C])
forces[0:-1]-=FEA(cumuMasses[0:-1], masses[1:C], radii[0:-1], radii[1,C])
forces[0]-=FFB(N[0], 0, radii[0]) #Special case because radii does not include 0
forces[1:C]-=FFA(N[1:], radii[0:-1], radii[1:C])
#print(forces)
I suspect that the culprit may be the variables I declared at the beginning, which are all probably Python's normal objects, but I don't know how to fix that.
The error is complaining about the "return arrays"; it doesn't say anything about the type of the "input arrays".
Working from @motto's comment, and previous SO with this error, I realized that this error has to do with the out parameter of a ufunc.
For example with a common np.add, if I call it with:
In [92]: np.add(1,2,3)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[92], line 1
----> 1 np.add(1,2,3)
np.add = <ufunc 'add'>
np = <module 'numpy' from 'C:\\Users\\paul\\miniconda3\\lib\\site-packages\\numpy\\__init__.py'>
TypeError: return arrays must be of ArrayType
The help signature of np.add is
add(x1, x2, /, out=None, *, ...
It can be called with 3 arguments (plus kwargs), but the third is the out, which must be None or an array; 3 is wrong.
Specify an out of the right size array, and it works:
In [96]: np.add(1,2,np.array(5))
Out[96]: array(3)
Or to define a frompyfunc with a wrong number of inputs as you do:
In [97]: def foo(a,b,c): # 3 arguments
...: print(a,b,c)
...: return a+b+c
...:
In [98]: f = np.frompyfunc(foo,2,1) # mistakenly using 2
In [99]: f(1,2,3) # and calling f with 3
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[99], line 1
----> 1 f(1,2,3)
f = <ufunc 'foo (vectorized)'>
TypeError: return arrays must be of ArrayType
Calling f with just 2 arguments gets around this out problem, but then runs into a problem with the number of arguments the foo requires:
In [103]: f(1,2)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[103], line 1
----> 1 f(1,2)
f = <ufunc 'foo (vectorized)'>
TypeError: foo() missing 1 required positional argument: 'c'
Specifying an array out, has another problem; frompyfunc returns an object dtype array:
In [104]: f(1,2,np.array(5))
---------------------------------------------------------------------------
UFuncTypeError Traceback (most recent call last)
Cell In[104], line 1
----> 1 f(1,2,np.array(5))
f = <ufunc 'foo (vectorized)'>
np = <module 'numpy' from 'C:\\Users\\paul\\miniconda3\\lib\\site-packages\\numpy\\__init__.py'>
UFuncTypeError: Cannot cast ufunc 'foo (vectorized)' output from dtype('O') to dtype('int32') with casting rule 'same_kind'
Using an object dtype array as out, gets us back to the problem providing enough arguments to foo:
In [105]: f(1,2,np.array(5,object))
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[105], line 1
----> 1 f(1,2,np.array(5,object))
f = <ufunc 'foo (vectorized)'>
np = <module 'numpy' from 'C:\\Users\\paul\\miniconda3\\lib\\site-packages\\numpy\\__init__.py'>
TypeError: foo() missing 1 required positional argument: 'c'
You are using frompyfunc because you need the r>x if test, which only works for scalars. Often there are better ways of dealing with tests like this.
Using frompyfunc:
In [112]: def foo(a,b):
...: if a>b:
...: return a+b
...: else:
...: return 0
In [113]: f = np.frompyfunc(foo,2,1)
In [115]: f(np.arange(5),2)
Out[115]: array([0, 0, 0, 5, 6], dtype=object)
Many ufunc like np.add handle a where/out pair of arguments:
In [116]: def foo1(a,b):
...: return np.add(a,b,where=a>b, out=np.zeros_like(a))
In [117]: foo1(np.arange(5),2)
Out[117]: array([0, 0, 0, 5, 6])
Or just the regular where:
In [119]: a=np.arange(5); b=2; np.where(a>b,a+b,0)
Out[119]: array([0, 0, 0, 5, 6])
This evaluates a+b in full, and then selects values; the [116] where is useful where the action produces an error or runtimewarning. You can also use masks to set some values according to condition, etc.
In short np.vectorze and np.frompyfunc can be useful when the function must take scalars. But even then a simple list comprehension might be just as good.
This old SO involving np.log could almost be used as a [duplicate]
TypeError: return arrays must be of ArrayType for a function that uses only floats
https://numpy.org/doc/stable/reference/generated/numpy.frompyfunc.html https://numpy.org/doc/stable/reference/ufuncs.html