Something went wrong. Try again.
[READ-ONLY] Mirror of https://github.com/emilyhunt/masters_project. My Master's project code, using a mixture density network to infer photometric redshifts to galaxies with spectroscopic redshifts as training data.
Something went wrong. Try again.
masters_project plotter_density.py
6.2 kB · 156 lines
Python
at master
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157"""Makes an SFR vs MStar density plot for the final report."""
import numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom mpl_toolkits.axes_grid1 import make_axes_locatable
# Make TeX labels work on plots#plt.rc('font', **{'family': 'serif', 'serif': ['DejaVu Sans']})#plt.rc('text', usetex=True)
# Constants and shitgrid_resolution = 20n_levels = 50n_shape = 60x_lim = [-10, 20]y_lim = [-10, 20]z_lim = [-2000, 50000]
# DATA READINGvalidation_no_spec_z = pd.read_csv('./final_run_data/validation_results_no_spec_z.csv')data_no_spec_z = pd.read_csv('./final_run_data/data_no_spec_z.csv')data_train = pd.read_csv('./final_run_data/data_validation.csv')
x = data_no_spec_z['sed_log_mass'].valuesy = data_no_spec_z['sed_log_sfr'].valuesz = np.asarray(np.abs(validation_no_spec_z['map'] - data_no_spec_z['z_phot_lit']) / (1 + 0.5*(validation_no_spec_z['map'] + data_no_spec_z['z_phot_lit'])))
# SET LIMITS OF TEST DATA# Set resolutionx_resolution = y_resolution = grid_resolution
# Setup the limits, ignoring anything not within said limitsallowed_x = np.where(np.logical_and(x_lim[0] < x, x < x_lim[1]), True, False)allowed_y = np.where(np.logical_and(y_lim[0] < y, y < y_lim[1]), True, False)allowed_z = np.where(np.logical_and(z_lim[0] < z, z < z_lim[1]), True, False)
# Now, only keep values that satisfy all conditionsallowed_values = np.logical_and(allowed_x, np.logical_and(allowed_y, allowed_z))x = x[allowed_values]y = y[allowed_values]z = z[allowed_values]
# TRAINING DATA LIMIT FINDING# Find the limits of the datax_train = data_train['sed_log_mass'].valuesy_train = data_train['sed_log_sfr'].values
allowed_x = np.where(np.logical_and(x_lim[0] < x_train, x_train < x_lim[1]), True, False)allowed_y = np.where(np.logical_and(y_lim[0] < y_train, y_train < y_lim[1]), True, False)allowed_values = np.logical_and(allowed_x, allowed_y)
x_train = x_train[allowed_values]y_train = y_train[allowed_values]
# Make an arbitrary shape of radius 1, then convert from polar coords back to cartesian# interior_angle = (n_shape - 2) * np.pi / n_shapeangles = np.linspace(0, 2 * np.pi - (2 * np.pi / n_shape), num=n_shape)shape_x = np.cos(angles)shape_y = np.sin(angles)
# Normalise the shape coords to be in the range [0, 1]shape_x = shape_x / (shape_x.max() - shape_x.min())shape_x -= shape_x.min()shape_y = shape_y / (shape_y.max() - shape_y.min())shape_y -= shape_y.min()
# Stretch the shape to fit across the datashape_x = (x_train.max() - x_train.min()) * shape_x + x_train.min()shape_y = (y_train.max() - y_train.min()) * shape_y + y_train.min()
# Make a big array of distances from each point of each point, then find the nearest one each timedistance_from_points = np.sqrt((shape_x.reshape(-1, 1) - x_train)**2 + (shape_y.reshape(-1, 1) - y_train)**2)point_ids = np.argmin(distance_from_points, axis=1)#point_ids = np.unique(point_ids)
# Grab the corresponding points to this from x and y train, and add the end on too so we plot a complete shapex_train_limits = x_train[point_ids]y_train_limits = y_train[point_ids]x_train_limits = np.append(x_train_limits, x_train_limits[0])y_train_limits = np.append(y_train_limits, y_train_limits[0])
# PLOTTING SHIT - TEST DATA# Calculate some ranges to grid overx_range = np.linspace(x.min(), x.max(), x_resolution)x_spacing = np.abs(x_range[1] - x_range[0])y_range = np.linspace(y.min(), y.max(), y_resolution)y_spacing = np.abs(y_range[1] - y_range[0])
# Make grid points for laterx_grid, y_grid = np.meshgrid(x_range, y_range)
# Cycle over the different grid points calculating the mean redshift in each placemean_z = np.zeros(x_grid.shape)training_data_density = np.zeros(x_grid.shape)i = 0for a_x in x_range: j = 0 for a_y in y_range: good_x = np.logical_and(a_x - 0.5*x_spacing < x, x < a_x + 0.5*x_spacing) good_y = np.logical_and(a_y - 0.5*y_spacing < y, y < a_y + 0.5*y_spacing) good_z = z[np.logical_and(good_x, good_y)]
good_x_train = np.logical_and(a_x - 0.5 * x_spacing < x_train, x_train < a_x + 0.5 * x_spacing) good_y_train = np.logical_and(a_y - 0.5 * y_spacing < y_train, y_train < a_y + 0.5 * y_spacing) training_data_density[i, j] = np.count_nonzero(np.logical_and(good_x_train, good_y_train))
# Only calculate the mean if there's more than one in this bin! if good_z.size > 0: mean_z[i, j] = np.mean(good_z) else: mean_z[i, j] = np.nan j += 1 i += 1
# PLOTTING TIME# The test datafig = plt.figure(figsize=(3, 3))ax = fig.add_subplot(1, 1, 1)#density_plot_object = ax.contourf(x_grid, y_grid, mean_z.T, n_levels, cmap=plt.get_cmap('plasma'), corner_mask=False)density_plot_object = ax.imshow(np.flip(mean_z.T, axis=0), extent=[x.min() - 0.5*x_spacing, x.max() + 0.5*x_spacing, y.min() - 0.5*y_spacing, y.max() + 0.5*y_spacing], cmap=plt.get_cmap('coolwarm'), alpha=1)
# Now, plot the training data's limitsax.contour(x_grid, y_grid, np.log(training_data_density.T + 1), 5, colors='k', linewidths=[2, 1, 1, 1, 1])#ax.plot(x_train_limits, y_train_limits, 'k-', label='Training data extent')#ax.plot(shape_x, shape_y, 'b-', label='Catching shape')#ax.plot(x_train, y_train, 'k.', alpha=0.1, ms=2)
# And let's label what the contours are#ax.text(1, 3, 'Training data\nnumber density', fontsize=8,# bbox=dict(boxstyle='round', ec=(0.0, 0.0, 0.0), fc=(1., 1.0, 1.0),))#ax.arrow(6, 1.5, 1, -0.5, head_width=0.3, head_length=0.3)
# Setup of the colorbardivider = make_axes_locatable(ax)cax = divider.append_axes('right', size='5%', pad=0.05)colorbar = fig.colorbar(density_plot_object, cax=cax, orientation='vertical')colorbar.set_label(r'$(z_{EAZY} - z_{ML}) \: / \: (1 + z_{EAZY})$')
# Labels and legendax.set_xlabel(r'$\log (M_{stars})$' + ' (stellar mass)')ax.set_ylabel(r'$\log (SFR)$' + ' (star formation rate)')#ax.legend(edgecolor='k', facecolor='w', fancybox=True, fontsize=8)
#fig.tight_layout()fig.savefig('./final_plots/pres_2_sfr_vs_mass.png', dpi=600, bbox_inches='tight')fig.show()