I'm running Python 2.7.11 from the Anaconda distribution (version 2.4.1)
on a MacBook Pro running Mac OS X version 10.11.2 (El Capitan)
I'm attempting to use numpy.ma.polyfit to perform a linear least square
fit on some data I have. I'm running NumPy version 1.10.1. I've observed
that in executing either numpy.polyfit or numpy.ma.polyfit I get the
following traceback:
/Users/user/anaconda/lib/python2.7/site-packages/numpy/lib/polynomial.py:594:
RankWarning: Polyfit may be poorly conditioned
warnings.warn(msg, RankWarning)
Traceback (most recent call last):
File "ComputeEnergy.py", line 132, in <module>
coeffs, covar = np.ma.polyfit( xfit, yfit, fit_degree,
rcond=rcondv, cov=True )
File
"/Users/user/anaconda/lib/python2.7/site-packages/numpy/ma/extras.py",
line 1951, in polyfit
return np.polyfit(x, y, deg, rcond, full, w, cov)
File
"/Users/user/anaconda/lib/python2.7/site-packages/numpy/lib/polynomial.py",
line 607, in polyfit
return c, Vbase * fac
ValueError: operands could not be broadcast together with shapes (6,6) (0,)
I've attached a stripped down version of the Python program I'm running.
Any suggestions?
Sam Dupree.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 07 22:25:29 2015
@author: user
"""
from astropy.time import Time
from datetime import datetime
import matplotlib.pyplot as plt
import numpy as np
import os
import string
# set physical constants
# -----------------------
# math constants
# --------------
rad_per_deg = np.pi / 180.
deg_per_rad = 1. / rad_per_deg
# orbital parameters for Earth
# ----------------------------
GMearth = 3.986004418e14 # in m**3/sec**2
Rearth = 6371.0e3 # in m
count = 0
x = []
y = []
fname = "TLE-ISS-25544.txt"
fhand = open(fname)
print
print "Processing TLEs in ", fname
print "------------------ "
print
for line in fhand:
count = count + 1
if line[0] == '1':
sat_number = line[2:7]
sat_classif = line[7:8]
intl_desig = line[9:17]
yr = int(line[18:20])
dayno = float(line[20:32])
if yr < 50:
yr = yr + 2000
else:
yr = yr + 1900
elif line[0] == '2':
mean_motion_kozai = float(line[52:63])
mean_motion = mean_motion_kozai * 2.e0 * np.pi / 86400.e0
total_eng = -0.5e0 * ( mean_motion * GMearth ) ** (2.e0/3.e0)
print
print ' sat number = ', sat_number, ' sat classif = ', sat_classif, ' Internatl Desig = ', intl_desig
print ' yr, dayno = ', yr, dayno
print ' mean motion (Revs/day) = ', mean_motion_kozai, ' mean motion (rad/sec) = ', mean_motion
print ' Total Energy (m^2/sec^2) = ', total_eng
print
x.append(dayno)
y.append(total_eng)
fhand.close()
count = len(x)
fit_span = 9
fit_cycles = count - fit_span + 1
fit_degree = 5
fit_offset = 0
evalpt = fit_span - 1
midpt = int( np.ceil(fit_span / 2) )
eps = 2.e-16
print
print " Fit Control Parameters "
print " ---------------------- "
print " count = ", count
print " eps = ", eps
print " evalpt = ", evalpt
print " midpt = ", midpt
print " fit_span = ", fit_span
print " fir_degree = ", fit_degree
print " fit_offset = ", fit_offset
print " fit_cycles = ", fit_cycles
print
xplot = []
yplot = []
yplot_resid = []
dyplot = []
xdyplot = []
for ifit in range( 0, fit_cycles ):
xfit = []
yfit = []
coeffs = []
resids = []
rank = []
sv = []
rrcond = []
covar = []
for i in range( 0, fit_span ):
k = i + fit_offset
xfit.append( x[k] )
yfit.append( y[k] )
rcondv = eps * len(xfit)
# coeffs = np.polyfit( xfit, yfit, fit_degree, rcond=rcondv )
# coeffs, resids, rank, sv, rrcond = np.polyfit( xfit, yfit, fit_degree, rcond=rcondv, full=True )
# coeffs, covar = np.polyfit( xfit, yfit, fit_degree, rcond=rcondv, cov=True )
coeffs, covar = np.ma.polyfit( xfit, yfit, fit_degree, rcond=rcondv, cov=True )
rows = len(covar)
xpoly = xfit[evalpt]
ypoly = np.polyval( coeffs, xfit[evalpt] )
ypoly_resid = ypoly - yfit[evalpt]
xplot.append( xpoly )
yplot.append( ypoly )
yplot_resid.append( ypoly_resid )
p = np.poly1d( coeffs )
dp = np.polyder( p )
dypoly = dp( xfit[midpt] )
xdyplot.append( xfit[midpt] )
dyplot.append( dypoly )
print
print "ifit = ", ifit
print " fit_offset, fit_offset+fit_size = ", fit_offset, fit_offset+fit_span
print
print " rcondv = ", rcondv
print
print " coeffs = ", coeffs
print
print " resids = ", resids
print
print " rank = ", rank
print
print " sv = ", sv
print
print " rrcond = ", rrcond
print
print " covar = "
print(covar)
print
print('\n'.join([' '.join(['{:4}'.format(item) for item in row])
for row in covar]))
print
print " xfit[evalpt], yfit[eval] = ", xfit[evalpt], yfit[evalpt]
print " xpoly , ypoly = ", xpoly, ypoly
print " xpoly , ypoly_resid = ", xpoly, ypoly_resid
print
print " xfit[midpt] , dypoly = ", xfit[midpt], dypoly
fit_offset = fit_offset + 1
# plot the total energy vs. modified Julian date
# ----------------------------------------------
plt.figure( 'Figure 1. Total Energy versus Time' )
plt.plot( x, y )
plt.title('Total Energy versus Time' )
plt.xlabel( 'Day Number -- 2015' )
plt.ylabel( 'Energy [meters^2 / second^2]' )
plt.savefig("Total Energy vs.Time.jpg")
plt.show()
# plot the fitted total energy vs. modified Julian date
# -----------------------------------------------------
plt.figure( 'Figure 2. Fitted Total Energy versus Time' )
plt.plot( xplot, yplot )
plt.title('Fitted Total Energy versus Time' )
plt.xlabel( 'Day Number -- 2015' )
plt.ylabel( 'Energy [meters^2 / second^2]' )
plt.savefig("Fitted Total Energy vs.Time.jpg")
plt.show()
# plot the fitted total energy vs. modified Julian date
# -----------------------------------------------------
plt.figure( 'Figure 3. Total Energy Residuals versus Time' )
plt.plot( xplot, yplot_resid )
plt.title('Total Energy Residual versus Time' )
plt.xlabel( 'Day Number -- 2015' )
plt.ylabel( 'Energy [meters^2 / second^2]' )
plt.savefig("Total Energy Residual vs.Time.jpg")
plt.show()
# plot the fitted total energy vs. modified Julian date
# -----------------------------------------------------
plt.figure( 'Figure 4. Fitted Total Energy Derivative versus Time' )
plt.plot( xdyplot, dyplot )
plt.title('Fitted Total Energy Derivative versus Time' )
plt.xlabel( 'Day Number -- 2015' )
plt.ylabel( 'Energy [meters^2 / second^3]' )
plt.savefig("Fitted Total Energy Derivative vs.Time.jpg")
plt.show()
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion
