In [19]:
%pylab inline
import numpy
import astrohack_projections
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 16, 16
plot_kwargs = {"interpolation": "nearest",
"cmap": "afmhot",
"origin": "lower"}
file_path = "./illustris_galaxies/cutout_83.hdf5"
illustris_gal = astrohack_projections.illustris_model_and_image(file_path)
illustris_gal.set_image_shape((30, 40))
for i in xrange(16):
xi_hat, eta_hat = astrohack_projections.choose_random_projection()
alpha, beta, gamma = numpy.random.uniform(0.0, 360.0, 3)
intensity = 10
scale = 0.2 * numpy.exp(numpy.random.uniform())
xshift = numpy.random.uniform(13., 16.)
yshift = numpy.random.uniform(18., 21.)
psf_size = 1
bg = 0.
kwargs = {'alpha':alpha, 'beta':beta, 'gamma':gamma, 'intensity':intensity, 'scale':scale, 'xshift': xshift, 'yshift': yshift, 'bg':0.0, 'psf_size':psf_size}
illustris_gal.set_image_parameters(**kwargs)
illustris_gal.construct_image(xi_hat, eta_hat)
plt.subplot(4, 4, i+1)
a, mid = numpy.percentile(illustris_gal.get_image(), [15., 85.])
ran = mid - a
vmin = mid - 5. * ran
vmax = mid + 20. * ran
plt.imshow(numpy.sqrt(illustris_gal.get_image()), **plot_kwargs)
plt.colorbar()
In [50]:
# create an album of 32 such objects with random projection + normaly distributed noise
# plot parameters
plt.rcParams['figure.figsize'] = 16, 16
plot_kwargs = {"interpolation": "nearest",
"cmap": "afmhot",
"origin": "lower"}
# initial galaxy data
file_path = "./illustris_galaxies/cutout_83.hdf5"
illustris_gal = astrohack_projections.illustris_model_and_image(file_path)
illustris_gal.set_image_shape((30, 40))
# album and PSF initialisation
album = astrohack_projections.album_and_model()
psf = astrohack_projections.mixture_of_gaussians(2)
psf.add_gaussian(1., numpy.array([0., 0.]), numpy.eye(2)*1.)
# galaxy model I fit with
basevar = 0.5 * numpy.eye(3)
gal_model = astrohack_projections.galaxy_model_3d()
v = [3., 0., 0.]
gal_model.add_gaussian(1.0, numpy.array([0., -1., 0.]), basevar + numpy.outer(v,v))
v = [-1., 3., 0.]
gal_model.add_gaussian(1.0, numpy.array([2., 1., 0.]), basevar + numpy.outer(v,v))
v = [1., 3., 0.]
gal_model.add_gaussian(1.0, numpy.array([-2., 1., 0.]), basevar + numpy.outer(v,v))
for i in xrange(32):
# projection parameters
xi_hat, eta_hat = astrohack_projections.choose_random_projection()
alpha, beta, gamma = numpy.random.uniform(0.0, 360.0, 3)
intensity = 150
scale = 0.18 * numpy.exp(numpy.random.uniform())
xshift = numpy.random.uniform(13., 16.)
yshift = numpy.random.uniform(18., 21.)
psf_size = 1
bg = 0.
# illustris galaxy
kwargs = {'alpha':alpha, 'beta':beta, 'gamma':gamma, 'intensity':intensity, 'scale':scale, 'xshift': xshift, 'yshift': yshift, 'bg':0.0, 'psf_size':psf_size}
illustris_gal.set_image_parameters(**kwargs)
illustris_gal.construct_image(xi_hat, eta_hat)
# image
image = astrohack_projections.image_and_model()
image.set_shape((30, 40))
image.set_psf(psf)
kwargs = {'alpha':alpha, 'beta':beta, 'gamma':gamma, 'intensity':intensity, 'scale':scale, 'xshift': xshift, 'yshift': yshift, 'bg':0.0}
image.set_parameters(**kwargs)
image.set_galaxy(gal_model)
image.set_ivar(numpy.ones(image.get_shape()))
image.set_data(illustris_gal.get_image() + numpy.random.normal(size=image.get_shape())/ numpy.sqrt(image.ivar))
# album
album.add_image(image)
print len(album)
In [56]:
def showme(album):
plt.rcParams['figure.figsize'] = 20, 80
plot_kwargs = {"interpolation": "nearest",
"cmap": "afmhot",
"origin": "lower"}
for i in xrange(len(album)):
image = album.get_all_images()[i]
plt.subplot(16, 4, 2*i+1)
vmin = -5. / numpy.sqrt(numpy.median(image.get_ivar())) # assumes bg = 0
vmax = -2. * vmin # assumes bg = 0
plt.imshow(numpy.nan_to_num(numpy.sqrt(image.get_data())), **plot_kwargs)
plt.colorbar()
plt.subplot(16, 4, 2*i+2)
plt.imshow(numpy.nan_to_num(numpy.sqrt(image.get_synthetic())), **plot_kwargs)
plt.colorbar()
showme(album)
In [57]:
# INITIALIZATION BLOCK!
import scipy.optimize as op
galpar0 = numpy.array([1., 2., 0., 0., 1., 1., 1., 0., 0., 0.,
1., 0., 2., 0., 1., 1., 1., 0., 0., 0.,
1., 0., 0., 2., 1., 1., 1., 0., 0., 0.])
imgpar0 = (10., 30., 30., 0.1, 0.5, 14., 20., 0.)
print "album before", album(galpar0)
for image in album:
print "image before", image(imgpar0)
showme(album)
plt.savefig("./fitting_pngs/illustris_083/initialisation.pdf")
In [58]:
# run the thing!
num_runs = 8
for i in xrange(num_runs):
# album optimisation
galpar0 = album.galaxy.get_parameters_vector()
print "album before", album(galpar0)
result = op.minimize(album, galpar0, method="Powell")
galpar = result['x']
print "album after", album(galpar)
showme(album)
plt.savefig("./fitting_pngs/illustris_083/run_%s_album.pdf" % str(2*i).zfill(2))
plt.close()
# image optimisation
for image in album:
imgpar0 = image.get_parameters_vector()
print "image before", image(imgpar0)
result = op.minimize(image, imgpar0)
imgpar = result['x']
print "image after", image(imgpar)
showme(album)
plt.savefig("./fitting_pngs/illustris_083/run_%s_image.pdf" % str(2*i+1).zfill(2))
plt.close()