-4

Why num_r1(x) and num_r2(x) type numpy.ndarray, but num_r(t) is type float? How can I keep num_r(t) type as array?

def num_r(t):
    for x in t:
        if x>tx:
            return num_r2(x)
        else:
            return num_r1(x)

Thank you!

The complete example is below

# -*- coding: utf-8 -*
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
from pylab import *

#### physical parameters
c = 2.998*10**10
hp = 6.626*10**-27
hb = 1.055*10**-27
kb = 1.381*10**-16
g = 6.673*10**-8
me = 9.109*10**-28
mp = 1.673*10**-24
q = 4.803*10**-10  #gausi
sigT = 6.652*10**-25

# The evolution of the characteristic frequencies
p = 2.5
E52 = 1
epsB_r = 1
epse_r = 1

D28 = 1 
n1 = 1.0
nu15 = 1*10**(-5)
r014 = 1
g42 = 1
delt12 =1
g4 = g42*10**2.5

E0 = E52*10**52
eta = g4
N0 = E0/(g4*mp*c**2)



p_tx = 3**(1./3)*2**(4./3)*mp**(-1./3)*c**(-5./3)
tx = p_tx*n1**(-1./3)*eta**(-8./3)

p_num_r1 = 2**(11./2)*7**(-2)*mp**(5./2)*me**(-3)*pi**(-1./2)*q*p_tx**(-6)*2**30*3**18*10**12

p_nuc_r1 = 2**(-33./2)*3**(-4)*10**(-4)*me*mp**(-3./2)*c**(-2)*sigT**(-2)*pi**(-1./2)*q

p_Fmax_r1 = 2**(15./2)*3**(9./2)*10**30*p_tx**(-3./2)*10**(-56)*me*mp**(1./2)*c**3*sigT*q**(-1)*2**(1./2)*3**(-1)

p_num_r2 = 2**(11./2)*7**(-2)*mp**(5./2)*me**(-3)*pi**(-1./2)*q*p_tx**(54./35)*(2**5*3**3*10**2)**(-54./35)

p_nuc_r2 = 2**(-13./2)*3**2*pi**(-1./2)*me*mp**(-3./2)*c**(-2)*sigT**(-2)*q*p_tx**(-74./35)*(2**5*3**3*10**2)**(4./35)

p_Fmax_r2 = 2**(1./2)*3**(-1)*pi**(-1./2)*me*mp**(1./2)*c**3*sigT*q**(-1)*10**(-56)



num_r1 = lambda t : p_num_r1*eta**18*((p-2)/(p-1))**2*epse_r**2*epsB_r**(1./2)*n1**(5./2)*t**6*E52**(-2)

nuc_r1 = lambda t : p_nuc_r1*eta**(-4)*epsB**(-3./2)*n1**(-3./2)*t**(-2)

Fmax_r1 = lambda t : p_Fmax_r1*N0**t**(3./2)*n1*eta**6*E52**(-1./2)*D28**(-2)*epsB_r**(1./2)

num_r2 = lambda t : p_num_r2*((p-2)/(p-1))**2*n1**(-74./35)*n1**(74./105)*eta**(592./105)*E52**(-74./105)

nuc_r2 = lambda t : p_nuc_r2*eta**(172./105)*t**(4./35)*n1**(-167./210)*epsB_r**(-3./2)

Fmax_r2 = lambda t : N0*eta**(62./105)*n1**(37./210)*epsB_r**(1./2)*t**(-34./35)*D28**(-2)

def fspe(t,u):
  if num_r(t)<nuc_r(t):
    return np.where(u<num_r(t),(u/num_r(t))**(1./3)*Fmax_r(t),np.where(u<nuc_r(t),(u/num_r(t))**(-(p-1.)/2)*Fmax_r(t),(u/nuc_r(t))**(-p/2)*(nuc_r(t)/num_r(t))**(-(p-1.)/2)*Fmax_r(t)))
  else:
    return np.where(u<nuc_r(t),(u/nuc_r(t))**(1./3)*Fmax_r(t),np.where(u<num_r(t),(u/nuc_r(t))**(-1./2)*Fmax_r(t),(u/num_r(t))**(-p/2)*(num_r(t)/nuc_r(t))**(-1.2)*Fmax_r(t)))


def num_r(t):
 for x in t:
  if x>tx:
    return num_r2(x)
  else:
    return num_r1(x)
def nuc_r(t):
 for x in t:
  if t>tx:
    return nuc_r2(x)
  else:
    return nuc_r1(x)
def Fmax_r(t):
 for x in t:
  if t>tx:
    return Fmax_r2(x)
  else:
    return Fmax_r1(x)


i= np.arange(-4,6,0.1)
t = 10**i
dnum   = [math.log10(mmm) for mmm in num_r(t)]
dnuc   = [math.log10(j) for j in nuc_r(t)]
nu_obs = [math.log(2.4*10**17,10) for a in i]
plt.figure('God Bless: Observable Limit')
plt.title(r'$\nu_{obs}$ and $\nu_c$ and $\nu_m$''\nComparation')
plt.xlabel('Time: log t')
plt.ylabel(r'log $\nu$')
plt.axvline(math.log10(tx))
plt.plot(i,nu_obs,'.',label=r'$\nu_{obs}$')
plt.plot(i,dnum,'D',label=r'$\nu_m$')
plt.plot(i,dnuc,'s',label=r'$\nu_c$')
plt.legend()
plt.grid(True)
plt.savefig("nu_obs.eps", dpi=120,bbox_inches='tight')
plt.show()

But thereś a Error

TypeError  Traceback (most recent call last)
<ipython-input-250-c008d4ed7571> in <module>()
     95     i= np.arange(-4,6,0.1)
     96     t = 10**i
--->     97     dnum   = [math.log10(mmm) for mmm in num_r(t)]

TypeError: 'float' object is not iterable

2 Answers2

1

You should write your function as:

def num_r_(x):
  if x > tx:
      return num_r2(x)
  else:
      return num_r1(x)

And then pass it through np.vectorize to lift it from float to float to np.array to np.array

num_r = np.vectorize(num_r_)

From Efficient evaluation of a function at every cell of a NumPy array

And then when you use it in:

dnum = [math.log10(mmm) for mmm in num_r(t)]

You should rather do:

dnum = np.log10(num_r(t))

That is to say don't use the functions from the math module. Use those from the np module as they can take np.array as well as float.

As:

i = np.arange(-4,6,0.1)
t = 10**i

results in t being a np.array

Community
  • 1
  • 1
Dan D.
  • 73,243
  • 15
  • 104
  • 123
1

So i is an array (arange); so is t (a math expression of i).

def num_r(t):
 for x in t:
  if x>tx:
    return num_r2(x)
  else:
    return num_r1(x)

You iterate on t. x is an element of t. You test it and pass it through num_r2 or num_r1, and return immediately. So only the 1st element t is being processed. Thus the error - num_r returns one value, not an array.

You need to write num_r in a way that processes all the values of t, not just the first. A simple, crude way is

def num_r(t):
    result = []
    for x in t:
       if x>tx:
          value = num_r2(x)
       else:
          value = num_r1(x)
       result.append(value)
    # result = np.array(result)
    return result

Now num_r should return a list the same length as t, and can be use in the list comprehension

[math.log10(mmm) for mmm in num_r(t)]

num_r could be written as a list comprehension:

[(num_r2(x) if x>tx else num_r1(x)) for x in t]

You could have it return an array instead of a list, but as long as you are using it in the list comprehension, there's no need. A list is just fine.

If it did return an array, then you could replace the list comprehension with a numpy log operation, e.g.

np.log10(num_r(t))

If num_r1 and num_r2 are written so they can take an array (looks off hand like they are, but I haven't tested), you could write

def num_r(t):
     ind = t>tx
     result = np.zeros_like(t)
     result[ind] = num_r2(t[ind])
     result[~ind] = num_r1(t[~ind])
     return result

The idea is to find a mask of a values in t that are >tx, and pass all those through num_r2 at once; similarly for num_r1; and collect the values in the correct slots of result. The result is an array that can be passed to np.log10. This should be quite a bit faster than iterating on t, or using np.vectorize.

There may be some errors in my suggestions, since I did not test them in an script or interpreter. But the underlying ideas should be correct and set you on the right path.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • While not actually faster in this case, `np.where(t>tx, num_r2(t), num_r1(t))` is a more concise way to express this. – Eric Mar 14 '16 at 05:19
  • Yes `where` should work. But practice with basic masked indexing is also desirable. – hpaulj Mar 14 '16 at 05:44
  • I really appreciate it ! Thanks for your help, though some other Errors occur, I believe your method is useful. I'd like to study the detail from your illustrations. Thank you again! God Bless You! – John Paul Qiang Chen Mar 14 '16 at 06:42