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






    

In [30]:

    

# initialisation
galpar0 = album_new.galaxy.get_parameters_vector()
print len(galpar0)
imgpar0 = (10., 30., 30., 0.1, 0.5, 14., 20., 0.)

print "album before", album_new(galpar0)
for image in album_new:
    print "image before", image(imgpar0)

showme(album_new)
plt.savefig("./fitting_pngs/illustris_242959_complex/initialisation.pdf")


    

    




100
album before 122395323.424
image before 7753875.73337
image before 6124179.15108
image before 7688199.81114
image before 784916.969444
image before 1172269.0661
image before 8728463.44021
image before 13999323.3387
image before 1392297.68871
image before 1007281.92019
image before 1979290.586
image before 2442506.24125
image before 3242410.23289
image before 1288348.29693
image before 2244332.66231
image before 6102391.52241
image before 3159549.64073
image before 1331344.52038
image before 915621.33846
image before 1517871.38907
image before 3663172.73136
image before 2188927.67713
image before 4113243.4138
image before 1971160.35588
image before 2581194.99315
image before 1907400.3404
image before 674153.436075
image before 6744445.16814
image before 10838710.7892
image before 1191658.22706
image before 1045946.71339
image before 7353816.60464
image before 5256770.59796






    




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






    

In [ ]:

    

# run the thing!
import scipy.optimize as op
num_runs = 20

for i in xrange(num_runs):
    # album optimisation
    galpar0 = album_new.galaxy.get_parameters_vector()
    print "album before", album_new(galpar0)
    result = op.minimize(album_new, galpar0, method="Powell")
    galpar = result['x']
    print "album after", album_new(galpar)
    
    showme(album_new)
    plt.savefig("./fitting_pngs/illustris_242959_complex/run_%s_album.pdf" % str(2*i).zfill(2))
    plt.close()
    
    # image optimisation
    for image in album_new:
        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_new)
    plt.savefig("./fitting_pngs/illustris_242959_complex/run_%s_image.pdf" % str(2*i+1).zfill(2))
    plt.close()


    

    




album before 





    




/Users/dalyabaron/.virtualenvs/astro/lib/python2.7/site-packages/scipy/optimize/optimize.py:1778: RuntimeWarning: invalid value encountered in double_scalars
  tmp2 = (x - v) * (fx - fw)
/Users/dalyabaron/.virtualenvs/astro/lib/python2.7/site-packages/scipy/optimize/optimize.py:1779: RuntimeWarning: invalid value encountered in double_scalars
  p = (x - v) * tmp2 - (x - w) * tmp1