%pylab inline
from scipy import stats
from scipy.optimize import curve_fit
tau = 2.0 * math.pi
with open('data/hygfull.csv') as f:
print f.read(250) + '...'
import csv
with open('data/hygfull.csv') as f:
stars = list(csv.DictReader(f))
print 'The database holds', len(stars), 'stars'
print 'An example star:', stars[123]
from spectral_classification import build_color_chart
color_chart = build_color_chart('starcolors.txt')
print 'The RGB color for class G2(V) is:\n',
print color_chart['G2(V)']
usefuls = [
star for star in stars
if star['Spectrum'] in color_chart and
star['ColorIndex'] != ''
]
print 'There are', len(usefuls), '/', len(stars), 'useful stars'
print 'Which is', 100.0 * float(len(usefuls)) / len(stars), '%'
color_index = np.array(
[ float(star['ColorIndex']) for star in usefuls ])
spectral_code = np.array(
[ star['Spectrum'] for star in usefuls ])
print 'First few spectral codes:', spectral_code[:5]
print 'First few color indexes:', color_index[:5]
rgb_color = np.array(
[ color_chart[code] for code in spectral_code ])
print 'First few RGB colors:\n', rgb_color[:5]
red = rgb_color[:,0]
green = rgb_color[:,1]
blue = rgb_color[:,2]
big = (8.0, 16.0)
scatter(color_index, red, s=1.0)
scatter(color_index, green, s=1.0)
scatter(color_index, blue, s=1.0)
result = stats.linregress(color_index, red)
slope, intercept, r_value, p_value, std_err = result
print slope, intercept, r_value, p_value, std_err
scatter(color_index, red, s=1.0)
x0, x1 = -1, 1.5
plot([x0, x1], [intercept + x0 * slope,
intercept + x1 * slope])
def fancy_regression(color):
result = stats.linregress(color_index[color < 1.0],
color[color < 1.0])
slope, intercept, r_value, p_value, std_err = result
print slope, intercept, r_value, p_value, std_err
scatter(color_index, color, s=1.0)
x0, x1 = axis()[:2]
plot([x0, x1], [intercept + x0 * slope,
intercept + x1 * slope])
axis([None, None, 0.3, 1.1])
fancy_regression(red)
fancy_regression(blue)
def func(x, a, b, c):
return a * exp(-b*x) + c
x = np.linspace(0,4,50)
y = func(x, 2.5, 1.3, 0.5)
yn = y + 0.2*np.random.normal(size=len(x))
popt, pcov = curve_fit(func, x, yn)
print popt
print pcov
def parabola(x, a, b, c):
return a*x*x + b*x + c
popt, pcov = curve_fit(parabola, color_index, green)
print popt
print pcov
scatter(color_index, green, s=1.0)
x0, x1, y0, y1 = axis()
xspace = np.linspace(x0, x1, 100)
yvalue = parabola(xspace, popt[0], popt[1], popt[2])
plot(xspace, yvalue)
axis([None, None, 0.6, 1.1])
def indexed_color(color_index):
x = color_index
r = 0.448834253080 * x + 0.697688812105
g = -0.22906088 * x * x + 0.38220694 * x + 0.77957636
b = -0.349557269948 * x + 1.170004977010
return np.array((r, g, b)).clip(0.0, 1.0)
print indexed_color(0.0)
print indexed_color(1.0)
print indexed_color(2.0)
xspace = np.linspace(-1.0, 3.0, 21)[2:-2]
colors = indexed_color(xspace).transpose()
scatter(xspace, xspace * 0.0 + 0.5,
s=100.0, c=colors)