In [1]:
%matplotlib inline
%config IPython.matplotlib.backend = 'retina'
%config InlineBackend.figure_format = 'retina'
import numpy as np
from astropy.table import Table, hstack, join
import dustmaps.bayestar, dustmaps.sfd
import astropy.units as units
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import tensorflow as tf
%cd ..
from starlight.models_cy import *
%cd Starlight
In [2]:
import gaia_tools.select
import gaia_tools.load as gload
In [3]:
tgas_cat = gload.tgas()
twomass = gload.twomass()
In [4]:
tgas_cat.dtype.names
Out[4]:
In [5]:
twomass.dtype.names
Out[5]:
In [6]:
tsf = gaia_tools.select.tgasSelect()
In [7]:
indx = tsf.determine_statistical(tgas_cat,
twomass['j_mag'],
twomass['k_mag'])
In [8]:
#import healpy
#healpy.mollview(title="")
#healpy.projplot(tgas_cat['l'][indx],
# tgas_cat['b'][indx],
# 'k,', lonlat=True, alpha=0.03);
In [9]:
#indx[:] = True
indx &= twomass['matched']
indx &= np.isfinite(twomass['j_msigcom'])
indx &= np.isfinite(twomass['h_msigcom'])
indx &= np.isfinite(twomass['k_msigcom'])
indx &= np.isfinite(tgas_cat['phot_g_mean_flux_error'])
indx &= np.isfinite(tgas_cat['parallax_error'])
indx &= tgas_cat['parallax_error'] > 0
indx &= (tgas_cat['parallax'] / tgas_cat['parallax_error'] > 1/1) # Main cut
print('Number of objects=', indx.sum())
In [10]:
bayestar = dustmaps.bayestar.BayestarQuery(max_samples=1)
sfd = dustmaps.sfd.SFDQuery()
In [403]:
tgas_cat.dtype
Out[403]:
In [422]:
In [323]:
sel = np.random.choice(np.where(indx)[0], int(1e6), replace=False)
sel = np.random.choice(np.where(indx)[0], indx.sum(), replace=False)
In [429]:
varpis = np.asarray(tgas_cat['parallax'][sel]).astype('f8')
varpis_var = np.asarray(tgas_cat['parallax_error'][sel]).astype('f8') ** 2.0
nobj = varpis.size
ls = np.asarray(tgas_cat['l'][sel]).astype('f8')
bs = np.asarray(tgas_cat['b'][sel]).astype('f8')
distances = (1000/varpis)
coords = SkyCoord(ls*units.deg, bs*units.deg, distance=distances*units.pc, frame='galactic')
ras, decs = coords.icrs.ra.rad, coords.icrs.dec.rad
ebv = bayestar(coords, mode='median')
ebv2 = sfd(coords)
ind2 = ~np.isfinite(ebv)
ebv[ind2] = 0#ebv2[ind2]
# Coefs from Lauren Anderson.
# Second set from http://astroquery.readthedocs.io/en/latest/irsa/irsa_dust.html
# see also http://iopscience.iop.org/article/10.1088/0004-637X/737/2/103/pdf
# for others bandpasses.
G_RedCoeff = 2.55
J_RedCoeff = 0.709 # 0.723
H_RedCoeff = 0.449 # 0.46
K_RedCoeff = 0.30 # 0.31
zeropoint = 25.52477
gmag = zeropoint - 2.5*np.log(np.asarray(tgas_cat['phot_g_mean_flux'][sel]).astype('f8'))/np.log(10.)
gmag_err = -2.5/np.log(10.) * np.asarray(tgas_cat['phot_g_mean_flux_error'][sel] / tgas_cat['phot_g_mean_flux'][sel]).astype('f8')
obsmags = np.asarray(twomass['j_mag'][sel]).astype('f8')
obscolors = np.vstack((
np.asarray(twomass['j_mag'][sel] - twomass['k_mag'][sel]).astype('f8'),
np.asarray(twomass['j_mag'][sel] - twomass['h_mag'][sel]).astype('f8'),
np.asarray(twomass['k_mag'][sel] - twomass['h_mag'][sel]).astype('f8')
)).T
#obscolors = np.asarray(twomass['j_mag'][sel] - twomass['k_mag'][sel]).astype('f8')[:, None]
ncols = obscolors.shape[1]
obsmagsandcolors = np.zeros((nobj, ncols+1))
obsmagsandcolors_var = np.zeros((nobj, ncols+1, ncols+1))
obsmagsandcolors[:, 0] = np.asarray(twomass['j_mag'][sel]).astype('f8')
obsmagsandcolors[:, 1] = np.asarray(twomass['j_mag'][sel] - twomass['k_mag'][sel]).astype('f8')
obsmagsandcolors[:, 2] = np.asarray(twomass['j_mag'][sel] - twomass['h_mag'][sel]).astype('f8')
obsmagsandcolors[:, 3] = np.asarray(twomass['k_mag'][sel] - twomass['h_mag'][sel]).astype('f8')
obsmagsandcolors_var[:, 0, 0] = np.asarray(twomass['j_msigcom'][sel]).astype('f8') ** 2.
if True:
obsmagsandcolors_var[:, 0, 1] = np.asarray(twomass['j_msigcom'][sel]).astype('f8') ** 2.
obsmagsandcolors_var[:, 0, 2] = np.asarray(twomass['j_msigcom'][sel]).astype('f8') ** 2.
obsmagsandcolors_var[:, 1, 2] = np.asarray(twomass['j_msigcom'][sel]).astype('f8') ** 2.
obsmagsandcolors_var[:, 2, 1] = np.asarray(twomass['j_msigcom'][sel]).astype('f8') ** 2.
obsmagsandcolors_var[:, 1, 0] = np.asarray(twomass['j_msigcom'][sel]).astype('f8') ** 2.
obsmagsandcolors_var[:, 2, 0] = np.asarray(twomass['j_msigcom'][sel]).astype('f8') ** 2.
obsmagsandcolors_var[:, 1, 1] = np.asarray(twomass['j_msigcom'][sel]**2. + twomass['k_msigcom'][sel]**2.).astype('f8')
obsmagsandcolors_var[:, 2, 2] = np.asarray(twomass['j_msigcom'][sel]**2. + twomass['h_msigcom'][sel]**2.).astype('f8')
obsmagsandcolors_var[:, 3, 3] = np.asarray(twomass['k_msigcom'][sel]**2. + twomass['h_msigcom'][sel]**2.).astype('f8')
dustcoefs = np.array([J_RedCoeff,
J_RedCoeff - K_RedCoeff,
J_RedCoeff - H_RedCoeff,
K_RedCoeff - H_RedCoeff])[:ncols+1]
magname = 'J'
colnames = ['J-K', 'J-H', 'K-H'][:ncols]
obsabsmags = obsmags + 5*np.log10(varpis) - 10
dustamps = 1*ebv.astype('f8')
varpimagvar = (varpis_var**0.5 * 5 / varpis / np.log(10))**2.0
In [430]:
fig, axs = plt.subplots(2, ncols, figsize=(10, 6), sharex=False, sharey=False)
if ncols == 1:
axs = axs.reshape((-1, 1))
rr = [[-0.5, 1.35], [-4, 6]]
nbins = 150
for i in range(ncols):
axs[0, i].hist2d(obsmagsandcolors[:, i+1],
obsmagsandcolors[:, 0] + 5*np.log10(varpis) - 10,
nbins, range=rr, cmap="gray_r")
axs[1, i].hist2d(obscolors[:, i] - dustamps * dustcoefs[i+1],
obsabsmags - dustamps * dustcoefs[0],
nbins, range=rr, cmap="gray_r")
for j in range(2):
axs[j, i].set_ylabel(magname)
axs[j, i].set_xlabel(colnames[i])
axs[j, i].set_ylim(rr[1][::-1])
axs[j, i].set_xlim(rr[0])
fig.tight_layout()
In [428]:
subsel = np.random.choice(nobj, int(2e2))
fig, axs = plt.subplots(2, ncols, figsize=(10, 6), sharex=False, sharey=False)
if ncols == 1:
axs = axs.reshape((-1, 1))
rr = [[-0.5, 1.35], [-4, 6]]
for i in range(ncols):
axs[0, i].errorbar(obsmagsandcolors[subsel, i+1],
obsmagsandcolors[subsel, 0] + 5*np.log10(varpis[subsel]) - 10,
xerr=obsmagsandcolors_var[subsel, i+1, i+1]**0.5,
yerr=(obsmagsandcolors_var[subsel, 0, 0] + varpimagvar[subsel])**0.5,
markersize=2, alpha=0.25, fmt="o", lw=1)
axs[1, i].errorbar(obsmagsandcolors[subsel, i+1] - dustamps[subsel] * dustcoefs[i+1],
obsmagsandcolors[subsel, 0] + 5*np.log10(varpis[subsel]) - 10 - dustamps[subsel] * dustcoefs[0],
xerr=obsmagsandcolors_var[subsel, i+1, i+1]**0.5,
yerr=(obsmagsandcolors_var[subsel, 0, 0] + varpimagvar[subsel])**0.5,
markersize=2, alpha=0.25, fmt="o", lw=1)
for j in range(2):
axs[j, i].set_ylabel(magname)
axs[j, i].set_xlabel(colnames[i])
axs[j, i].set_ylim(rr[1][::-1])
axs[j, i].set_xlim(rr[0])
fig.tight_layout()
In [327]:
import scipy.linalg
num = 10000
rhos = np.zeros((num, 3))
rhos_out = np.zeros((num, 3))
for i in range(num):
rhos[i, :] = np.random.uniform(0, 2e-0, 3)
#skewsym = np.array([[0, rhos[i, 0]], [-rhos[i, 0], 0]])
skewsym = np.zeros((3, 3))
skewsym[np.tril_indices(3, k=-1)] = rhos[i, :]
skewsym[np.triu_indices(3, k=1)] = -rhos[i, :]
ortho = scipy.linalg.expm(skewsym)
#ortho2 = np.diag(np.repeat(1, 3)) + skewsym + np.dot(skewsym, skewsym) / 2 + np.dot(np.dot(skewsym, skewsym), skewsym) / 6
#print(ortho)
#print(ortho2)
vvars = np.random.uniform(0, 1, 3)
#vvars[1] = vvars[0]*2
covar = np.dot(ortho.T, np.dot(np.diag(vvars), ortho))
rhos_out[i, 0] = covar[1, 0] / np.sqrt(covar[0, 0]*covar[1, 1])
rhos_out[i, 1] = covar[2, 0] / np.sqrt(covar[0, 0]*covar[2, 2])
rhos_out[i, 2] = covar[2, 1] / np.sqrt(covar[1, 1]*covar[2, 2])
fig, axs = plt.subplots(1, 3, figsize=(10, 3))
axs[0].hist2d(rhos_out[:, 0], rhos_out[:, 1], cmap="Greys");
axs[1].hist2d(rhos_out[:, 1], rhos_out[:, 2], cmap="Greys");
axs[2].hist2d(rhos_out[:, 0], rhos_out[:, 2], cmap="Greys");
In [328]:
sigma = 3
nmarggrid = 5
varpimins = varpis - sigma * varpis_var**0.5
varpimins[varpimins<1e-4] = 1e-4
varpimaxs = varpis + sigma * varpis_var**0.5
deltavarpis = (varpimaxs - varpimins) / (nmarggrid-1)
varpigrids = varpimins[:, None] + deltavarpis[:, None] * np.arange(nmarggrid)[None, :]
logvarpigridsp10_padded = np.zeros((nobj, nmarggrid, ncols + 1))
logvarpigridsp10_padded[:, :, 0] = 5*np.log10(varpigrids) - 10.
In [337]:
obsabsmagcolors_minmaxs = np.zeros((ncols+1, 2))
obsabsmagcolors_minmaxs[0, 0] = np.min(obsabsmags)
obsabsmagcolors_minmaxs[1:, 0] = np.min(obscolors, axis=0)
obsabsmagcolors_minmaxs[0, 1] = np.max(obsabsmags)
obsabsmagcolors_minmaxs[1:, 1] = np.max(obscolors, axis=0)
ngridpoints = 15
lines = [np.linspace(obsabsmagcolors_minmaxs[i, 0], obsabsmagcolors_minmaxs[i, 1], ngridpoints) for i in range(ncols+1)]
mids = [(line[1:] + line[:-1])/2 for line in lines]
flatmeshes = [mesh.ravel() for mesh in np.meshgrid(*mids, indexing="ij")]
H, E = np.histogramdd(np.hstack((obsabsmags[:, None], obscolors)), bins=lines)
selection = np.where(H > 1)
In [338]:
fig, axs = plt.subplots(1, ncols, figsize=(10, 3), sharex=False, sharey=False)
if ncols == 1:
axs = np.array([axs])
for i in range(ncols):
ind = tuple([j for j in range(ncols+1) if j != 0 and j != i+1])
post = H.sum(axis=ind)
post = np.where(post > nobj//1000, 1, 0)
axs[i].pcolormesh(E[i+1], E[0], post, cmap="Greys")
ind = np.random.choice(nobj, 100, replace=False)
axs[i].scatter(obscolors[ind, i], obsabsmags[ind], s=1, color='red')
axs[i].set_ylim([obsabsmagcolors_minmaxs[0, 1], obsabsmagcolors_minmaxs[0, 0]])
axs[i].set_ylabel(magname)
axs[i].set_xlabel(colnames[i])
In [339]:
nbins = 40
binamps = np.repeat(1./nbins, nbins)
if False: # select positions among objects
sel2 = np.random.choice(nobj, nbins, replace=False)
binmus = np.hstack((obsabsmags[sel2][:, None], obscolors[sel2, :]))
else: # select positions uniformly on grid
selection = np.where(H.ravel() > nobj // 100)[0]
sel2 = np.random.choice(selection, nbins, replace=True)
binmus = np.vstack([mesh[sel2] for mesh in flatmeshes]).T
binvars = np.vstack((np.random.uniform(0.1, 1.0, nbins)**2,
[np.random.uniform(0.03, 0.1, nbins)**2 for i in range(ncols)])).T
print(binmus)
bincovars = np.zeros((nbins, ncols+1, ncols+1))
for b in range(nbins):
bincovars[b, :, :] = np.diag(binvars[b, :])
from matplotlib.patches import Ellipse
fig, axs = plt.subplots(1, ncols, figsize=(10, 3), sharex=False, sharey=False)
if ncols == 1:
axs = np.array([axs])
rr = [[-0.5, 1.35], [-4, 6]]
for i in range(ncols):
axs[i].hist2d(obscolors[:, i] - dustamps * dustcoefs[i+1],
obsabsmags - dustamps * dustcoefs[0],
150, range=rr, cmap="gray_r")
for b in range(nbins):
bincovar = np.matrix([[bincovars[b, 0, 0], bincovars[b, i+1, 0]],
[bincovars[b, 0, i+1], bincovars[b, i+1, i+1]]])
w, v = np.linalg.eig(bincovar[:, :])
angle = 90. + 180/np.pi * 0.5 * np.arctan(2 * bincovar[0, 1] / (bincovar[1, 1] - bincovar[0, 0]))
e = Ellipse(xy=binmus[b, :][np.array([i+1, 0])],
width=4*w[0]**0.5, height=4*w[1]**0.5, angle=angle, ec='k', fill=False)
axs[i].add_artist(e)
e.set_alpha(binamps[b]**0.5)
axs[i].scatter(binmus[b, i+1], binmus[b, 0], c='y')
axs[i].scatter(binmus[b, i+1], binmus[b, 0], c='y')
axs[i].set_ylabel(magname)
axs[i].set_xlabel(colnames[i])
axs[i].set_ylim(rr[1][::-1])
axs[i].set_xlim(rr[0])
fig.tight_layout()
In [370]:
def fun(zs_in):
zs = 1e-2 + 0.99 / (1 + np.exp(-zs_in))
fac = np.array([1.0]*nbins)
zsb = np.array([1.0]*nbins)
for i in range(nbins-1):
fac[i] = 1. - zs[i]
zsb[i+1] = zs[i]
binamps = np.cumprod(zsb) * fac
return binamps
def obj(zs_in):
binamps_target = np.repeat(1/nbins, nbins)
return np.sum((fun(zs_in) - binamps_target)**2.)
from scipy.optimize import minimize
res = minimize(obj, x0=np.random.uniform(-2, 2, nbins-1))
print(res.success, res.x)
binamps = fun(res.x)
print(binamps)
zs = res.x
In [371]:
import tensorflow as tf
from scipy.special import logit
def tfzs_to_alphas(zs):
fac = tf.concat((1 - zs, tf.ones((1, 1, 1))), axis=1)
zsb = tf.concat((tf.ones((1, 1, 1)), zs), axis=1)
fs = tf.cumprod(zsb, axis=1) * fac
return fs
Zs = tf.Variable(zs.reshape((1, nbins-1, 1)), dtype=tf.float32, name="Zs")
Binmus_in = tf.Variable(logit((binmus+4)/10), dtype=tf.float32, name="Binmus_in")
Binvars_in = tf.Variable(logit((binvars-1e-4)/1.5), dtype=tf.float32, name="Binvars_in")
Binmus = -4 + 10 * tf.sigmoid(Binmus_in)
Binvars = 1e-4 + 1.5 * tf.sigmoid(Binvars_in)
Binamps = tfzs_to_alphas(1e-2 + 0.99*tf.sigmoid(Zs))
Obsmagsandcolors = tf.placeholder(shape=(None, ncols + 1), dtype=tf.float32)
Obsmagsandcolors_var = tf.placeholder(shape=(None, ncols + 1, ncols + 1), dtype=tf.float32)
if True:
Dustamps = tf.placeholder(shape=(None, ), dtype=tf.float32)
else:
Dust_lmn_real = tf.Variable(dust_lmn[:, :, None, None].real, dtype=tf.float32, name="Dust_lmn_real")
Dust_lmn_imag = tf.Variable(dust_lmn[:, :, None, None].imag, dtype=tf.float32, name="Dust_lmn_imag")
Dust_lmn = tf.complex(Dust_lmn_real, Dust_lmn_imag)
SFB_mat = tf.placeholder(shape=(Lt, N, Npix, Nr), dtype=tf.complex64)
Dust_map = tf.cumsum(tf.exp(tf.real(tf.reduce_sum(SFB_mat * Dust_lmn, axis=(0, 1))) - 4.0), axis=1)
#Dust_map_diff = tf.Variable(np.random.randn(Npix, Nr), dtype=tf.float32, name="Dust_map_diff")
#Dust_map = tf.cumsum(tf.sigmoid(Dust_map_diff)/Nr, axis=1)
AngRadlocsgrid = tf.placeholder(shape=(None, nmarggrid, 2), dtype=tf.int32)
Dustamps = tf.gather_nd(Dust_map, AngRadlocsgrid)
Dustcoefs = tf.placeholder(shape=(ncols+1, ), dtype=tf.float32)
Varpigrids = tf.placeholder(shape=(None, nmarggrid, ), dtype=tf.float32)
Logvarpigridsp10_padded = tf.placeholder(shape=(None, nmarggrid, ncols + 1), dtype=tf.float32)
Deltavarpis = tf.placeholder(shape=(None, ), dtype=tf.float32)
Obsvarpis = tf.placeholder(shape=(None, ), dtype=tf.float32)
Obsvarpis_var = tf.placeholder(shape=(None, ), dtype=tf.float32)
if True:
Skewsym_vals = tf.Variable(1e-2+np.zeros(((nbins*(ncols+1)*ncols//2), )),
dtype=tf.float32, name="skewsym_vals")
inds = [np.vstack([np.repeat(b, (ncols+1)*ncols//2), np.tril_indices(ncols+1, k=-1)]) for b in range(nbins)]
indices = np.concatenate([list(zip(*inds[b])) for b in range(nbins)])
Indices = tf.constant([list(i) for i in indices], dtype=tf.int64)
#Skewsymt2 = tf.sigmoid(tf.sparse_to_dense(sparse_indices=Indices, output_shape=[nbins, ncols+1, ncols+1], \
# sparse_values=Skewsym_vals, default_value=0, \
# validate_indices=True))
Skewsym = 2 * tf.sigmoid(tf.scatter_nd(Indices, Skewsym_vals, [nbins, ncols+1, ncols+1]))
Skewsym -= tf.transpose(Skewsym, perm=[0, 2, 1])
fac = 1.0
if False:
Newterm = tf.concat([tf.diag(tf.ones((ncols+1)))[None, :, :] for b in range(nbins)], axis=0)
Orthogonal = 1*Newterm
for i in range(1, 10):
fac *= i
Newterm = tf.matmul(Newterm, Skewsym)
Orthogonal += Newterm / fac
else:
Orthogonal = tf.concat([tf.diag(tf.ones((ncols+1)))[None, :, :] for b in range(nbins)], axis=0)\
+ Skewsym + tf.matmul(Skewsym, Skewsym) / 2 + tf.matmul(tf.matmul(Skewsym, Skewsym), Skewsym) / 6
Binvarsdiag = tf.concat([tf.diag(Binvars[b, :])[None, :, :] for b in range(nbins)], axis=0)
Bincovars = tf.matmul(tf.transpose(Orthogonal, perm=[0, 2, 1]),
tf.matmul(Binvarsdiag, Orthogonal))
Bincovars = tf.scalar_mul(0.5, Bincovars + tf.transpose(Bincovars, perm=[0, 2, 1]))
else:
Skewsym_vals = tf.constant(np.random.randn((nbins*(ncols+1)*ncols//2), ),
dtype=tf.float32, name="skewsym_vals")
Bincovars = tf.concat([tf.diag(Binvars[b, :])[None, :, :] for b in range(nbins)], axis=0)
Covmatrices = Bincovars[None, :, :, :] + Obsmagsandcolors_var[:, None, :, :] # nobj, nbins, ncols+1, ncols+1
Covmatrices_chol = tf.cholesky(Covmatrices)[:, :, None, :, :] *\
tf.ones((1, nbins, nmarggrid, ncols+1, ncols+1)) # nobj, nbins, nmarggrid, ncols+1, ncols+1
Logdet_mags = tf.log(tf.matrix_determinant(Covmatrices))
Logdet_varpis = tf.log(Obsvarpis_var[:, None])
# nobj, nbins, nk, ncols + 1, 1
Meanvecs = Obsmagsandcolors[:, None, None, :, None] \
- Dustcoefs[None, None, None, :, None] * Dustamps[:, None, None, None, None] \
- Binmus[None, :, None, :, None] + Logvarpigridsp10_padded[:, None, :, :, None]
Chi2s_mags = tf.reduce_sum(
tf.multiply(Meanvecs, tf.cholesky_solve(Covmatrices_chol, Meanvecs)), axis=(3, 4)) # nobj, nbins
# for debugging:
#Covmatrices = Covmatrices[:, :, None, :, :] * tf.ones((1, nbins, nmarggrid, ncols+1, ncols+1))
#Chi2s_mags = tf.reduce_sum(
# tf.multiply(Meanvecs, tf.matrix_solve(Covmatrices, Meanvecs)), axis=(3, 4)) # nobj, nbins
Chi2s_varpis = tf.divide(tf.pow(Obsvarpis[:, None] - Varpigrids, 2), Obsvarpis_var[:, None])
# nobj, nbins, nk
Logprobgrid = tf.log(Deltavarpis[:, None, None] * Binamps) +\
tf.scalar_mul(-0.5, Chi2s_mags + Logdet_mags[:, :, None] + Chi2s_varpis[:, None, :] + Logdet_varpis[:, None, :])
MinusLnprob = - tf.reduce_sum(tf.reduce_logsumexp(Logprobgrid, axis=(1, 2)), axis=0)
In [372]:
Optimizer = tf.train.AdamOptimizer(learning_rate=1e-2)\
.minimize(MinusLnprob, var_list=[Zs, Binmus_in, Binvars_in, Skewsym_vals])#, Dust_lmn_real, Dust_lmn_imag]) #
# GradientDescentOptimizer AdadeltaOptimizer AdamOptimizer
with tf.Session() as sess:
sess.run(tf.global_variables_initializer())
for i in range(3001):
ind = np.random.choice(nobj, nobj // 500, replace=False)
_, mlnp, binamps2, binmus2, binvars2, bincovars2, skewsym_vals, skewsym =\
sess.run([Optimizer, MinusLnprob, Binamps, Binmus, Binvars, Bincovars, Skewsym_vals, Skewsym],
feed_dict={Obsmagsandcolors: obsmagsandcolors[ind, :],
Obsmagsandcolors_var: obsmagsandcolors_var[ind, :, :],
Dustamps: dustamps[ind], Dustcoefs: dustcoefs, Varpigrids: varpigrids[ind, :],
Logvarpigridsp10_padded: logvarpigridsp10_padded[ind, :, :],
Deltavarpis: deltavarpis[ind], Obsvarpis: varpis[ind], Obsvarpis_var: varpis_var[ind]
})
if i % 100 == 0:
print(i, mlnp, end=" ")
In [382]:
from matplotlib.patches import Ellipse
fig, axs = plt.subplots(3, ncols, figsize=(10, 10), sharex=False, sharey=False)
if ncols == 1:
axs = axs.reshape((-1, 1))
rr = [[-0.5, 1.35], [-4, 6]]
nbs2 = np.random.multinomial(nobj, binamps2.ravel()/(binamps2.sum()+1e-10))
cmap = "bone_r"
points = np.vstack([np.random.multivariate_normal(binmus2[b, :], bincovars2[b, :, :], size=nbs2[b])
for b in range(nbins)])
for i in range(ncols):
for k in range(2):
axs[k, i].hist2d(obscolors[:, i] - dustamps * dustcoefs[i+1],
obsabsmags - dustamps * dustcoefs[0],
100, range=rr, cmap=cmap)
axs[2, i].hist2d(points[:, i+1], points[:, 0], 100, range=rr, cmap=cmap)
for b in np.where(binamps2.ravel() > np.max(binamps2)/100)[0]:#range(nbins):
bincovar = np.matrix([[bincovars2[b, 0, 0], bincovars2[b, i+1, 0]],
[bincovars2[b, 0, i+1], bincovars2[b, i+1, i+1]]])
w, v = np.linalg.eig(bincovar[:, :])
angle = 90. + 180/np.pi * 0.5 * np.arctan(2 * bincovar[0, 1] / (bincovar[1, 1] - bincovar[0, 0]))
e = Ellipse(xy=binmus2[b, :][np.array([i+1, 0])], width=4*w[0]**0.5, height=4*w[1]**0.5, angle=angle, ec='k', fill=False)
axs[1, i].add_artist(e)
e.set_alpha(binamps2[0, b, 0]**0.5)
axs[1, i].scatter(binmus2[b, i+1], binmus2[b, 0], c='y')
for k in range(3):
axs[k, i].set_ylim(rr[1][::-1])
axs[k, i].set_xlim(rr[0])
axs[k, i].set_ylabel(magname)
axs[k, i].set_xlabel(colnames[i])
fig.tight_layout()
fig.savefig('multicolor_hrd_gmm.png', dpi=300)
In [ ]: