We want to use the Radio Galaxy Zoo click data as training data. To do this, we first need to convert the raw click data to a useful label — most likely the $(x, y)$ coordinate of the host galaxy for a given set of radio emissions.
Let's start by analysing what we're looking at in the raw data.
In [79]:
import collections
import operator
from pprint import pprint
import sqlite3
import sys
import warnings
import matplotlib.pyplot
import numpy
import scipy.stats
import sklearn.mixture
%matplotlib inline
sys.path.insert(1, '..')
import crowdastro.data
import crowdastro.show
warnings.simplefilter('ignore', UserWarning) # astropy always raises warnings on Windows.
We'll grab a random subject (complete, from the ATLAS survey) and look at its classifications.
In [4]:
subject = crowdastro.data.db.radio_subjects.find_one({'metadata.survey': 'atlas', 'state': 'complete',
'zooniverse_id': 'ARG0003r18'})
In [5]:
crowdastro.show.subject(subject)
matplotlib.pyplot.show()
This is a kinda nice subject to look at because there are multiple radio sources. This is often true for Radio Galaxy Zoo subjects. Let's look at a classification.
In [10]:
classification = crowdastro.data.db.radio_classifications.find_one({'subject_ids': subject['_id']})
pprint(classification)
We only really care about the annotations, so I've pulled those out. Some of them aren't useful — the langauge, user agent, and start/end times here, though there may potentially be more in other classifications. It looks like the ones that are useful contain both 'ir' and 'radio' keys. I'll use "annotations" to refer to these exclusively from now on.
Each annotation contains an IR component and a radio component. The IR component is obviously where people clicked, though the format is a little unusual — it seems to have the option for multiple clicks, though I can't come up with a situation in which this would arise. I think we can safely ignore all but the first IR coordinate. The radio component describes which contours are included in this annotation. Volunteers may click on multiple contours that they think have been emitted from the same radio source, and these are grouped here. The combination of radio contours here is important as it affects what the volunteers believe they are classifying (see Banfield et al. 2015 for more details).
We need to pull out something describing which combination of radio contours were selected, as well as the actual IR coordinate associated with each one.
In [7]:
def make_radio_combination_signature(radio_annotation):
"""Generates a unique signature for a radio annotation.
radio_annotation: 'radio' dictionary from a classification.
-> Something immutable
"""
# My choice of immutable object will be a tuple of the xmax values,
# sorted to ensure determinism, and rounded to nix floating point errors.
# Note that the x scale is not the same as the IR scale, but the scale factor is
# included in the annotation, so I have multiplied this out here for consistency.
# Sometimes, there's no scale, so I've included a default scale.
xmaxes = [round(float(c['xmax']) * float(c.get('scale_width', '2.1144278606965172')), 14)
for c in radio_annotation.values()]
return tuple(sorted(xmaxes))
def read_classification(classification):
"""Converts a raw RGZ classification into radio combinations and IR locations.
classification: RGZ classification dictionary.
-> dict mapping radio combination signatures to IR locations.
"""
result = {}
for annotation in classification['annotations']:
if 'radio' not in annotation:
# This is a metadata annotation and we can ignore it.
continue
radio_signature = make_radio_combination_signature(annotation['radio'])
if annotation['ir'] == 'No Sources':
ir_location = None
else:
ir_location = (float(annotation['ir']['0']['x']), float(annotation['ir']['0']['y']))
result[radio_signature] = ir_location
return result
In [8]:
read_classification(classification)
Out[8]:
Now let's try running that over a few classifications, to see if anything breaks.
In [9]:
for classification in crowdastro.data.db.radio_classifications.find().limit(100):
print(read_classification(classification))
This seems to work pretty well! The next step is to freeze the whole classification database in this simplified form. I'll do this in a script in crowdastro.__main__.
This section is summarised in crowdastro.labels.parse_classification.
Now we can move on to using these to find the consensus for each subject.
Given these simplified data points, it should be a lot easier to figure out what the consensus is for a given subject: Pull up all the classifications for a given subject, find the most common full radio signature, and get the clicks for each resulting radio signature.
I've frozen the previous results into a neat little database, so we should be good to go.
In [61]:
conn = sqlite3.connect('../crowdastro-data/processed.db')
conn.row_factory = sqlite3.Row
cur = conn.cursor()
classifications = list(cur.execute('SELECT full_radio_signature, part_radio_signature, source_x, source_y '
'FROM classifications WHERE subject_id=?', [str(subject['_id'])]))
frs_counter = collections.Counter([c['full_radio_signature'] for c in classifications])
most_common_frs = frs_counter.most_common(1)[0][0]
radio_consensus_classifications = collections.defaultdict(list)
for classification in classifications:
if classification['full_radio_signature'] == most_common_frs:
radio_consensus_classifications[classification['part_radio_signature']].append((classification['source_x'],
classification['source_y']))
conn.close()
In [62]:
radio_consensus_classifications
Out[62]:
Pulling out the coordinates for each radio source was reasonably straightforward. Now, we want to use these to figure out which location these volunteers agree on. Let's just focus on \xe3\x0c}? for now. Plotting it to get an idea of what we have:
In [76]:
crowdastro.show.ir(subject)
xs = [a[0] * crowdastro.config.get('click_to_fits_x') for a in radio_consensus_classifications[b'\xe3\x0c}?']]
ys = [crowdastro.config.get('fits_image_height') - a[1] * crowdastro.config.get('click_to_fits_y')
for a in radio_consensus_classifications[b'\xe3\x0c}?']]
matplotlib.pyplot.xlim((170, 200))
matplotlib.pyplot.ylim((170, 185))
matplotlib.pyplot.scatter(xs, ys, marker='+')
matplotlib.pyplot.show()
A Gaussian assumption is not reasonable — maybe the volunteers think there are two galaxies that the radio source could belong to? I think a Gaussian mixture assumption is reasonable, though, so I will use PG-means to try and find this.
In [146]:
def pg_means(points, significance=0.01, projections=24):
"""Cluster points with the PG-means algorithm."""
k = 1
while True:
# Fit a Gaussian mixture model with k components.
gmm = sklearn.mixture.GMM(n_components=k)
try:
gmm.fit(points)
except ValueError:
return None
for _ in range(projections):
# Project the data to one dimension.
projection_vector = numpy.random.random(size=(2,))
projected_points = points @ projection_vector
# Project the model to one dimension.
# We need the CDF in one dimension, so we'll sample some data points and project them.
n_samples = 1000
samples = gmm.sample(n_samples) @ projection_vector
samples.sort()
def cdf(x):
for sample, y in zip(samples, numpy.arange(n_samples) / n_samples):
if sample >= x:
break
return y
_, p_value = scipy.stats.kstest(projected_points, numpy.vectorize(cdf))
if p_value < significance:
# Reject the null hypothesis.
break
else:
# Null hypothesis was not broken.
return gmm
k += 1
Let's try that.
In [114]:
xs = [a[0] * crowdastro.config.get('click_to_fits_x') for a in radio_consensus_classifications[b'\xe3\x0c}?']]
ys = [crowdastro.config.get('fits_image_height') - a[1] * crowdastro.config.get('click_to_fits_y')
for a in radio_consensus_classifications[b'\xe3\x0c}?']]
points = numpy.vstack([numpy.array(xs), numpy.array(ys)])
crowdastro.show.ir(subject)
matplotlib.pyplot.xlim((170, 200))
matplotlib.pyplot.ylim((170, 185))
matplotlib.pyplot.scatter(xs, ys, marker='+')
matplotlib.pyplot.scatter(*pg_means(points.T, significance=0.02).means_.T)
matplotlib.pyplot.show()
That seems alright. Let's try with the other radio source.
In [145]:
xs = [a[0] * crowdastro.config.get('click_to_fits_x')
for a in radio_consensus_classifications[b'\xb4\xd7\x1c?']
if a[0] is not None]
ys = [crowdastro.config.get('fits_image_height') - a[1] * crowdastro.config.get('click_to_fits_y')
for a in radio_consensus_classifications[b'\xb4\xd7\x1c?']
if a[1] is not None]
points = numpy.vstack([numpy.array(xs), numpy.array(ys)])
crowdastro.show.ir(subject)
matplotlib.pyplot.xlim((80, 120))
matplotlib.pyplot.ylim((80, 150))
matplotlib.pyplot.scatter(xs, ys, marker='+')
matplotlib.pyplot.scatter(*pg_means(points.T, significance=0.02, projections=24).means_.T)
matplotlib.pyplot.show()
Note that we get different results for significance = 0.01 and significance = 0.02. This is unfortunately a hyperparameter that will need to be tweaked. Let's bundle this up into a function.
In [150]:
def get_subject_consensus(subject, conn=None, significance=0.02):
conn.row_factory = sqlite3.Row
cur = conn.cursor()
classifications = list(cur.execute('SELECT full_radio_signature, part_radio_signature, source_x, source_y '
'FROM classifications WHERE subject_id=?', [str(subject['_id'])]))
frs_counter = collections.Counter([c['full_radio_signature'] for c in classifications])
most_common_frs = frs_counter.most_common(1)[0][0]
radio_consensus_classifications = collections.defaultdict(list)
for classification in classifications:
if classification['full_radio_signature'] == most_common_frs:
radio_consensus_classifications[classification['part_radio_signature']].append((classification['source_x'],
classification['source_y']))
consensus = {} # Maps radio signatures to (x, y) NumPy arrays.
for radio_signature in radio_consensus_classifications:
xs = [a[0] * crowdastro.config.get('click_to_fits_x')
for a in radio_consensus_classifications[radio_signature]
if a[0] is not None]
ys = [crowdastro.config.get('fits_image_height') - a[1] * crowdastro.config.get('click_to_fits_y')
for a in radio_consensus_classifications[radio_signature]
if a[1] is not None]
points = numpy.vstack([xs, ys])
gmm = pg_means(points.T, significance=significance, projections=24)
consensus[radio_signature] = gmm.means_[gmm.weights_.argmax()]
return consensus
In [152]:
conn = sqlite3.connect('../crowdastro-data/processed.db')
consensus = get_subject_consensus(subject, conn=conn)
conn.close()
In [161]:
matplotlib.pyplot.figure(figsize=(10, 10))
matplotlib.pyplot.axis('off')
crowdastro.show.subject(subject)
matplotlib.pyplot.scatter(*numpy.array(list(consensus.values())).T)
matplotlib.pyplot.show()
Similarly, I'll run this over the whole dataset and freeze the results. The consensus algorithm is located in crowdastro.labels.find_subject_consensus and I have also added the freezing script to crowdastro.__main__.