[Scipy-tickets] [SciPy] #280: Failure of n=1 case in gen_roots_and_weights

SciPy scipy-tickets at scipy.net
Sun Oct 8 18:04:39 CDT 2006


#280: Failure of n=1 case in gen_roots_and_weights
---------------------------+------------------------------------------------
 Reporter:  jcarlson       |       Owner:  somebody
     Type:  defect         |      Status:  new     
 Priority:  normal         |   Milestone:          
Component:  scipy.special  |     Version:  devel   
 Severity:  normal         |    Keywords:          
---------------------------+------------------------------------------------
 Changes since 0.5.1 have broken the orthogonal polynomial functions in the
 n=1 case. For instance:

 {{{
 >>> from scipy.special import hermite
 >>> hermite(0)
 poly1d([ 1.])
 >>> hermite(1)
 Traceback (most recent call last):
   File "<stdin>", line 1, in <module>
   File "/home/jc/lib/python/scipy/special/orthogonal.py", line 316, in
 hermite
     p = orthopoly1d(x,w,hn,kn,wfunc,(-inf,inf),monic)
   File "/home/jc/lib/python/scipy/special/orthogonal.py", line 78, in
 __init__
     equiv_weights = [weights[k] / wfunc(roots[k]) for k in
 range(len(roots))]
 TypeError: object of type 'numpy.complex128' has no len()
 >>> hermite(2)
 poly1d([  4.00000000e+00,  -8.88178420e-16,  -2.00000000e+00])
 }}}


 The problem can be traced back to the function gen_roots_and_weights() in
 scipy.special.orthogonal (lines 111-113):

 {{{
     sortind = x.real.argsort()
     answer.append(x[sortind])
     answer.append((mu*v[0]**2)[sortind])
 }}}

 If {{{len(x) == 1}}}, then {{{sortind = array([0])}}}, and
 {{{x[sortind]}}} gives the ''value'' {{{x[0]}}} rather than the ''array''
 {{{x}}}.  This leads to trouble, as seen in the error output above.  A
 very naive patch is included below.

 {{{
 Index: Lib/special/orthogonal.py
 ===================================================================
 --- Lib/special/orthogonal.py   (revision 2247)
 +++ Lib/special/orthogonal.py   (working copy)
 @@ -108,9 +108,13 @@
                  np.diagflat(sqrt_bn,1) +
                  np.diagflat(sqrt_bn,-1)))
      answer = []
 -    sortind = x.real.argsort()
 -    answer.append(x[sortind])
 -    answer.append((mu*v[0]**2)[sortind])
 +    if n <= 1:
 +        answer.append(x)
 +        answer.append(mu*v[0]**2)
 +    else:
 +        sortind = x.real.argsort()
 +        answer.append(x[sortind])
 +        answer.append((mu*v[0]**2)[sortind])
      return answer

  # Jacobi Polynomials 1               P^(alpha,beta)_n(x)
 }}}

-- 
Ticket URL: <http://projects.scipy.org/scipy/scipy/ticket/280>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list