"TypeError: return arrays must be of ArrayType" even though all input arrays have the same type

Question

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.

Answer 1

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