In [43]:
%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 = "/Users/dalyabaron/Downloads/cutout_242959.hdf5"
illustris_gal = astrohack_projections.illustris_model_and_image(file_path)
illustris_gal.set_image_shape((60, 80))

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 = 1.5
    scale = 0.015 * numpy.exp(numpy.random.uniform())
    xshift = numpy.random.uniform(29., 31.)
    yshift = numpy.random.uniform(39., 41.)
    psf_size = 1.5
    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()
    
    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(illustris_gal.get_image()+1, **plot_kwargs)


Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['gamma', 'beta']
`%matplotlib` prevents importing * from pylab and numpy

In [41]:
# 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 = "/Users/dalyabaron/Downloads/cutout_242959.hdf5"
illustris_gal = astrohack_projections.illustris_model_and_image(file_path)
illustris_gal.set_image_shape((60, 80))

# 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 = numpy.random.uniform(0, 3, size=3) #[3., 0., 0.]
gal_model.add_gaussian(1.0, numpy.array([0., -1., 0.]), basevar + numpy.outer(v,v))
v = numpy.random.uniform(0, 3, size=3) #[-1., 3., 0.]
gal_model.add_gaussian(1.0, numpy.array([2., 1., 0.]), basevar + numpy.outer(v,v))
v = numpy.random.uniform(0, 3, size=3) #[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 = 20
    scale = 0.015 * numpy.exp(numpy.random.uniform())
    xshift = numpy.random.uniform(29.5, 31.)
    yshift = numpy.random.uniform(39.5, 41.)
    psf_size = 1.5
    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((60, 80))
    image.set_psf(psf)
    scale = 0.5
    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)


32

In [42]:
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)


/Users/dalyabaron/.virtualenvs/astro/lib/python2.7/site-packages/ipykernel/__main__.py:11: RuntimeWarning: invalid value encountered in sqrt

In [18]:
plt.rcParams['figure.figsize'] = 12, 8

# plot for example two of these galaxies and their light histogram
image_1 = album.get_all_images()[0]
image_2 = album.get_all_images()[1]

sum_1_x = numpy.sum(image_1.get_data(), axis=0)
sum_1_y = numpy.sum(image_1.get_data(), axis=1)
sum_2_x = numpy.sum(image_2.get_data(), axis=0)
sum_2_y = numpy.sum(image_2.get_data(), axis=1)

plt.subplot(2, 3, 1)
plt.imshow(numpy.nan_to_num(numpy.sqrt(image_1.get_data())), **plot_kwargs)
plt.subplot(2, 3, 2)
plt.plot(sum_1_x)
plt.xlim(0, 80)
plt.subplot(2, 3, 3)
plt.plot(sum_1_y)
plt.xlim(0, 80)

plt.subplot(2, 3, 4)
plt.imshow(numpy.nan_to_num(numpy.sqrt(image_2.get_data())), **plot_kwargs)
plt.subplot(2, 3, 5)
plt.plot(sum_2_x)
plt.xlim(0, 80)
plt.subplot(2, 3, 6)
plt.plot(sum_2_y)
plt.xlim(0, 80)


/Users/dalyabaron/.virtualenvs/astro/lib/python2.7/site-packages/ipykernel/__main__.py:13: RuntimeWarning: invalid value encountered in sqrt
/Users/dalyabaron/.virtualenvs/astro/lib/python2.7/site-packages/ipykernel/__main__.py:22: RuntimeWarning: invalid value encountered in sqrt
Out[18]:
(0, 80)

In [27]:
basevar = 0.1 * numpy.eye(3)
gal_model = astrohack_projections.galaxy_model_3d()
num_of_gaussians = 10

for i in xrange(num_of_gaussians):
    mu = numpy.random.uniform(10, 50, size=3)
    v = numpy.random.uniform(-3, 3, size=3)
    fi = basevar + numpy.outer(v, v)
    gal_model.add_gaussian(1.0, mu, fi)

# album and PSF initialisation
album_new = astrohack_projections.album_and_model()
psf = astrohack_projections.mixture_of_gaussians(2)
psf.add_gaussian(1., numpy.array([0., 0.]), numpy.eye(2)*1.)

for i in xrange(32):
    image = album.get_all_images()[i]
    data = image.get_data()
    parameters = image.get_parameters()
    print parameters
    # image
    image = astrohack_projections.image_and_model()
    image.set_shape((60, 80))
    image.set_psf(psf)
    parameters['scale'] = 1
    image.set_parameters(**parameters)
    image.set_galaxy(gal_model)
    image.set_ivar(numpy.ones(image.get_shape()))
    image.set_data(data)
    # album
    album_new.add_image(image)
    
showme(album_new)
album_new.set_galaxy(gal_model)
plt.savefig("./fitting_pngs/illustris_242959_complex/initialisation.pdf")


{'beta': 17.286919886785078, 'scale': 1, 'bg': 0.0, 'yshift': 40.948658007306726, 'xshift': 30.66280774934001, 'alpha': 176.92318705358673, 'intensity': 20, 'gamma': 339.68804315247417}
{'beta': 335.0037339994243, 'scale': 1, 'bg': 0.0, 'yshift': 39.63285243952303, 'xshift': 30.9735348724383, 'alpha': 7.9292522577579083, 'intensity': 20, 'gamma': 162.56933906924331}
{'beta': 179.08784303333007, 'scale': 1, 'bg': 0.0, 'yshift': 39.66586099626983, 'xshift': 29.9488218852435, 'alpha': 55.369875754731197, 'intensity': 20, 'gamma': 39.767653528878661}
{'beta': 137.82220813174257, 'scale': 1, 'bg': 0.0, 'yshift': 39.681374694522866, 'xshift': 30.13277480960852, 'alpha': 129.01729822354679, 'intensity': 20, 'gamma': 3.2425764038465488}
{'beta': 12.96039184721109, 'scale': 1, 'bg': 0.0, 'yshift': 39.62138366665347, 'xshift': 29.6652055251155, 'alpha': 215.59838912553147, 'intensity': 20, 'gamma': 249.71244588409795}
{'beta': 308.248051762996, 'scale': 1, 'bg': 0.0, 'yshift': 40.86505169460173, 'xshift': 30.13104106495362, 'alpha': 5.9402140026045425, 'intensity': 20, 'gamma': 24.811683207066793}
{'beta': 248.37148489898988, 'scale': 1, 'bg': 0.0, 'yshift': 39.55207831721144, 'xshift': 30.82677121363938, 'alpha': 188.72237792825641, 'intensity': 20, 'gamma': 63.849782053768529}
{'beta': 259.2689459942157, 'scale': 1, 'bg': 0.0, 'yshift': 40.51662423248716, 'xshift': 30.22283181494849, 'alpha': 21.408586733274845, 'intensity': 20, 'gamma': 249.28282427344953}
{'beta': 160.98692136624587, 'scale': 1, 'bg': 0.0, 'yshift': 40.89134067382771, 'xshift': 29.866576576614882, 'alpha': 314.93254919082193, 'intensity': 20, 'gamma': 164.10124075097215}
{'beta': 293.28950577580201, 'scale': 1, 'bg': 0.0, 'yshift': 40.403134437861695, 'xshift': 30.52672123720768, 'alpha': 0.2643000840886689, 'intensity': 20, 'gamma': 151.92268076911111}
{'beta': 297.59256294651448, 'scale': 1, 'bg': 0.0, 'yshift': 40.950439283028956, 'xshift': 29.950839379100323, 'alpha': 331.08552191458159, 'intensity': 20, 'gamma': 185.08318384451493}
{'beta': 242.91438782515755, 'scale': 1, 'bg': 0.0, 'yshift': 39.92057596410796, 'xshift': 30.692063309675504, 'alpha': 224.72727424414259, 'intensity': 20, 'gamma': 8.0090116833642266}
{'beta': 101.04914130920042, 'scale': 1, 'bg': 0.0, 'yshift': 40.68317681653718, 'xshift': 29.991035221743914, 'alpha': 24.129956525366755, 'intensity': 20, 'gamma': 120.96317247796152}
{'beta': 334.89433306918028, 'scale': 1, 'bg': 0.0, 'yshift': 39.57977796709863, 'xshift': 30.733199620043987, 'alpha': 187.08175272346506, 'intensity': 20, 'gamma': 138.57579235820108}
{'beta': 222.43903112850231, 'scale': 1, 'bg': 0.0, 'yshift': 39.8642367153819, 'xshift': 30.35202631783676, 'alpha': 357.29197968218551, 'intensity': 20, 'gamma': 206.89325526487724}
{'beta': 67.126923769603238, 'scale': 1, 'bg': 0.0, 'yshift': 40.34857512591552, 'xshift': 30.649030137352618, 'alpha': 124.86030645534326, 'intensity': 20, 'gamma': 301.28962339989812}
{'beta': 20.279760106786991, 'scale': 1, 'bg': 0.0, 'yshift': 40.01965594683611, 'xshift': 30.95212708821218, 'alpha': 122.05664414217617, 'intensity': 20, 'gamma': 16.075759611824555}
{'beta': 45.400287146098044, 'scale': 1, 'bg': 0.0, 'yshift': 40.86763118265568, 'xshift': 29.859546221687935, 'alpha': 328.4604600854716, 'intensity': 20, 'gamma': 232.19913532357828}
{'beta': 116.87572051751442, 'scale': 1, 'bg': 0.0, 'yshift': 40.40918485691773, 'xshift': 30.00120349889623, 'alpha': 100.72257216031394, 'intensity': 20, 'gamma': 354.73450214381194}
{'beta': 230.01233874226617, 'scale': 1, 'bg': 0.0, 'yshift': 40.07858887033961, 'xshift': 29.681670127190944, 'alpha': 327.10116142633149, 'intensity': 20, 'gamma': 239.78523455461402}
{'beta': 296.47652298916955, 'scale': 1, 'bg': 0.0, 'yshift': 39.897561384724234, 'xshift': 29.544578224716766, 'alpha': 24.571857588260464, 'intensity': 20, 'gamma': 27.011396759851031}
{'beta': 117.22130989336108, 'scale': 1, 'bg': 0.0, 'yshift': 40.78564921458812, 'xshift': 29.721879075335945, 'alpha': 151.36554494864478, 'intensity': 20, 'gamma': 59.418878822052704}
{'beta': 78.037847054638831, 'scale': 1, 'bg': 0.0, 'yshift': 40.05790772030172, 'xshift': 30.184627134513182, 'alpha': 169.95681182048685, 'intensity': 20, 'gamma': 183.56206044975718}
{'beta': 154.99676341299653, 'scale': 1, 'bg': 0.0, 'yshift': 40.4699893962817, 'xshift': 29.62603779087648, 'alpha': 238.92599778247398, 'intensity': 20, 'gamma': 183.57192264874578}
{'beta': 216.53429643618998, 'scale': 1, 'bg': 0.0, 'yshift': 40.82278988463932, 'xshift': 30.964435512525867, 'alpha': 293.79271725057083, 'intensity': 20, 'gamma': 19.479304070911802}
{'beta': 116.70937894031825, 'scale': 1, 'bg': 0.0, 'yshift': 39.63162658676126, 'xshift': 29.79576825110621, 'alpha': 45.416084291037876, 'intensity': 20, 'gamma': 247.2334552217612}
{'beta': 58.37617675459154, 'scale': 1, 'bg': 0.0, 'yshift': 39.634884930831625, 'xshift': 30.095486223537, 'alpha': 194.57420528218884, 'intensity': 20, 'gamma': 328.26686406825252}
{'beta': 180.89086222565822, 'scale': 1, 'bg': 0.0, 'yshift': 40.104672178826966, 'xshift': 30.586729075269265, 'alpha': 18.46176177392374, 'intensity': 20, 'gamma': 67.869179423044315}
{'beta': 221.68575389172176, 'scale': 1, 'bg': 0.0, 'yshift': 40.141979384746755, 'xshift': 30.486055905255594, 'alpha': 15.082766690369116, 'intensity': 20, 'gamma': 150.93133324018893}
{'beta': 343.96825118599816, 'scale': 1, 'bg': 0.0, 'yshift': 40.8297522910469, 'xshift': 30.60117514492223, 'alpha': 121.47010313409568, 'intensity': 20, 'gamma': 324.43427874735079}
{'beta': 22.434137968610223, 'scale': 1, 'bg': 0.0, 'yshift': 39.960281790503856, 'xshift': 30.68984162221386, 'alpha': 328.02378798505208, 'intensity': 20, 'gamma': 282.13865993814596}
{'beta': 259.7964124023888, 'scale': 1, 'bg': 0.0, 'yshift': 40.17830187652026, 'xshift': 29.76933337043577, 'alpha': 206.8338963467595, 'intensity': 20, 'gamma': 281.63376174073704}
/Users/dalyabaron/.virtualenvs/astro/lib/python2.7/site-packages/ipykernel/__main__.py:11: RuntimeWarning: invalid value encountered in sqrt