# standard first block for defining the token and makinhg it available as a system variable for the session
# token must be replaced with new one once it has expired
with open('/home/idies/keystone.token', 'r') as f:
token = f.read().rstrip('\n')# replace with your own token ID
import sys
sys.argv.append("--ident="+token)
import SciServer.CasJobs as CasJobs
from io import StringIO
import numpy as np
import pandas
from matplotlib import pyplot as plt
from sklearn.neighbors import KNeighborsRegressor
from astroML.plotting import scatter_contour
# some special settings
# ensure columns get written completely in notebook
pandas.set_option('display.max_colwidth', -1)
# do *not* show python warnings
import warnings
warnings.filterwarnings('ignore')
NOBJECTS = 20000
GAL_COLORS_DTYPE = [('u', float),
('g', float),
('r', float),
('i', float),
('z', float),
('redshift', float),
('redshift_err', float)]
Now we do this
# get data from CasJobs
query_text = ('\n'.join(
("SELECT TOP %i" % NOBJECTS,
" p.u, p.g, p.r, p.i, p.z, s.z, s.zerr",
"FROM PhotoObj AS p",
" JOIN SpecObj AS s ON s.bestobjid = p.objid",
"WHERE ",
" p.u BETWEEN 0 AND 19.6",
" AND p.g BETWEEN 0 AND 20",
" AND (s.class = 'GALAXY' OR s.class = 'QSO')")))
response = CasJobs.executeQuery(query_text, "DR12", token = token)
# read the result as a
output = StringIO(response.read().decode())
data = np.loadtxt(output, delimiter=',',skiprows=1, dtype=GAL_COLORS_DTYPE)
# show some of the data
print(data[1:10])
n_neighbors = 1
N = len(data)
# shuffle data
np.random.seed(0)
np.random.shuffle(data)
# put colors in a matrix
X = np.zeros((N, 4))
X[:, 0] = data['u'] - data['g']
X[:, 1] = data['g'] - data['r']
X[:, 2] = data['r'] - data['i']
X[:, 3] = data['i'] - data['z']
z = data['redshift']
# divide into training and testing data
Ntrain = N / 2
Xtrain = X[:Ntrain]
ztrain = z[:Ntrain]
Xtest = X[Ntrain:]
ztest = z[Ntrain:]
knn = KNeighborsRegressor(n_neighbors, weights='uniform')
zpred = knn.fit(Xtrain, ztrain).predict(Xtest)
rms = np.sqrt(np.mean((ztest - zpred) ** 2))
print("RMS error = %.2g" % rms)
axis_lim = np.array([-0.1, 2.5])
plt.figure(figsize=(12, 8))
ax = plt.axes()
plt.scatter(ztest, zpred, c='k', lw=0, s=4)
plt.plot(axis_lim, axis_lim, '--k')
plt.plot(axis_lim, axis_lim + rms, ':k')
plt.plot(axis_lim, axis_lim - rms, ':k')
plt.xlim(axis_lim)
plt.ylim(axis_lim)
plt.text(0.99, 0.02, "RMS error = %.2g" % rms,
ha='right', va='bottom', transform=ax.transAxes,
bbox=dict(ec='w', fc='w'), fontsize=16)
plt.title('Photo-z: Nearest Neigbor Regression')
plt.xlabel(r'$\mathrm{z_{spec}}$', fontsize=20)
plt.ylabel(r'$\mathrm{z_{phot}}$', fontsize=20)
plt.show()
Improve on this photo-z predictor.