<jupyter_start><jupyter_text>Fitting the Brick 23 IR CMD<jupyter_code>%matplotlib inline
%config InlineBackend.figure_format='retina'
# %config InlineBackend.figure_format='svg'

import os
import time
from glob import glob
import numpy as np

brick = 23
STARFISH = os.getenv("STARFISH")
isoc_dir = "b23ir_isoc"
lib_dir = "b23ir_lib"
synth_dir = "b23ir_synth"
fit_dir = "b23ir_fit"
wfc3_bands = ['F110W', 'F160W']<jupyter_output><empty_output><jupyter_text>Download IsochronesWe use a log-space age grid for ages less than a billion years, and a linear grid of every 0.5 Gyr thereafter.<jupyter_code>from astropy.coordinates import Distance
import astropy.units as u

from padova import AgeGridRequest, IsochroneRequest
from starfisher import LibraryBuilder

z_grid = [0.012, 0.015, 0.019, 0.024, 0.03]
delta_gyr = 0.5
late_ages = np.log10(np.arange(1e9 + delta_gyr, 13e9, delta_gyr * 1e9))
print late_ages
if not os.path.exists(os.path.join(STARFISH, isoc_dir)):
    for z in z_grid:
        # Young ages in log-grid
        r = AgeGridRequest(z,
                           min_log_age=6.6,
                           max_log_age=9.,
                           delta_log_age=0.5,
                           phot='wfc3', photsys_version='odfnew')
        for isoc in r.isochrone_set:
            isoc.export_for_starfish(os.path.join(STARFISH, isoc_dir),
                                     bands=wfc3_bands)
        
        # Old ages in linear grid
        for logage in late_ages:
            r = IsochroneRequest(z, logage,
                                 phot='wfc3', photsys_version='odfnew')
            r.isochrone.export_for_starfish(os.path.join(STARFISH, isoc_dir),
                                     bands=wfc3_bands)

d = Distance(785 * u.kpc)
builder = LibraryBuilder(isoc_dir, lib_dir,
                         nmag=len(wfc3_bands),
                         dmod=d.distmod.value,
                         iverb=3)
if not os.path.exists(builder.full_isofile_path):
    builder.install()<jupyter_output>[  9.           9.17609126   9.30103      9.39794001   9.47712125
   9.54406804   9.60205999   9.65321251   9.69897      9.74036269
   9.77815125   9.81291336   9.84509804   9.87506126   9.90308999
   9.92941893   9.95424251   9.97772361  10.          10.0211893
  10.04139269  10.06069784  10.07918125  10.09691001]<jupyter_text>Build the Isochrone Library and Synthesize CMD planes<jupyter_code>from collections import namedtuple
from starfisher import Lockfile
from starfisher import Synth
from starfisher import ExtinctionDistribution
from starfisher import ExtantCrowdingTable
from starfisher import ColorPlane
from m31hst.phatast import PhatAstTable

if not os.path.exists(os.path.join(STARFISH, synth_dir)):
    os.makedirs(os.path.join(STARFISH, synth_dir))

# No binning in our lockfile
lockfile = Lockfile(builder.read_isofile(), synth_dir, unbinned=True)

# No extinction, yet
young_av = ExtinctionDistribution()
old_av = ExtinctionDistribution()
rel_extinction = np.ones(len(wfc3_bands), dtype=float)
for av in (young_av, old_av):
    av.set_uniform(0.)

# Use PHAT AST from the outer field
crowd_path = os.path.join(synth_dir, "crowding.dat")
full_crowd_path = os.path.join(STARFISH, crowd_path)
tbl = PhatAstTable()
tbl.write_crowdfile_for_field(full_crowd_path, 0,
                              bands=('f110w', 'f160w'))
crowd = ExtantCrowdingTable(crowd_path)

# Define CMD planes
Lim = namedtuple('Lim', 'x y')
ir_lim = Lim(x=(-1., 2.), y=(25., 16.))
ir_cmd = ColorPlane((wfc3_bands.index('F110W'),
                     wfc3_bands.index('F160W')),
                    wfc3_bands.index('F160W'),
                    ir_lim.x,
                    (min(ir_lim.y), max(ir_lim.y)),
                    28.,
                    suffix='f110f160',
                    x_label=r'$\mathrm{F110W}-\mathrm{F160W}$',
                    y_label=r'$\mathrm{F110W}$',
                   dpix=0.1)
colour_planes = [ir_cmd]

synth = Synth(synth_dir, builder, lockfile, crowd,
              rel_extinction,
              young_extinction=young_av,
              old_extinction=old_av,
              planes=colour_planes,
              mass_span=(0.08, 150.),
              nstars=10000000)
if len(glob(os.path.join(STARFISH, synth_dir, "z*"))) == 0:
    synth.run_synth(n_cpu=4)
    synth.plot_all_hess(os.path.join(STARFISH, synth_dir, 'hess'))<jupyter_output>nx 30
ny 90
x.shape (90, 30)<jupyter_text>Export the dataset for StarFISH<jupyter_code>from astropy.table import Table
from m31hst import phat_v2_phot_path

if not os.path.exists(os.path.join(STARFISH, fit_dir)):
    os.makedirs(os.path.join(STARFISH, fit_dir))

data_root = os.path.join(fit_dir, "b23ir.")
full_data_path = os.path.join(STARFISH, '{0}f110f160'.format(data_root))
brick_table = Table.read(phat_v2_phot_path(brick), format='fits')

# Only use stars within the fitting box
c = brick_table['f110w_vega'] - brick_table['f160w_vega']
m = brick_table['f160w_vega']
sel = np.where((c > min(ir_lim.x)) & (c < max(ir_lim.x)) &
               (m > min(ir_lim.y)) & (m < max(ir_lim.y)))[0]
brick_table = brick_table[sel]

print("Fitting {0:d} stars".format(len(brick_table)))

if not os.path.exists(full_data_path):
    phot_dtype = np.dtype([('x', np.float), ('y', np.float)])
    photdata = np.empty(len(brick_table), dtype=phot_dtype)
    photdata['x'][:] = brick_table['f110w_vega'] - brick_table['f160w_vega']
    photdata['y'][:] = brick_table['f160w_vega']    
    np.savetxt(full_data_path, photdata, delimiter=' ', fmt='%.4f')
import matplotlib as mpl
import matplotlib.pyplot as plt
from androcmd.plot import contour_hess

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
contour_hess(ax, brick_table['f110w_vega'] - brick_table['f160w_vega'],
             brick_table['f160w_vega'], ir_lim.x, (max(ir_lim.y), min(ir_lim.y)),
             plot_args={'ms': 3})
ir_cmd.plot_mask(ax)
ax.set_xlabel(r'$\mathrm{F110W}-\mathrm{F160W}$')
ax.set_ylabel(r'$\mathrm{F110W}$')
ax.set_xlim(ir_lim.x)
ax.set_ylim(ir_lim.y)
fig.show()<jupyter_output>/Users/jsick/.virtualenvs/starfish/lib/python2.7/site-packages/matplotlib/contour.py:931: UserWarning: linewidths is ignored by contourf
  warnings.warn('linewidths is ignored by contourf')
/Users/jsick/.virtualenvs/starfish/lib/python2.7/site-packages/matplotlib/figure.py:387: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "<jupyter_text>Run StarFISH SFH<jupyter_code>from starfisher import SFH, Mask

mask = Mask(colour_planes)
sfh = SFH(data_root, synth, mask, fit_dir)
if not os.path.exists(sfh.full_outfile_path):
    sfh.run_sfh()
sfh_table = sfh.solution_table()
print(sfh_table)
from starfisher.sfhplot import SFHCirclePlot

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)

cp = SFHCirclePlot(sfh_table)
cp.plot_in_ax(ax, max_area=800)
for logage in np.log10(np.arange(1, 13, 1) * 1e9):
    ax.axvline(logage, c='0.8', zorder=-1)

ax.set_ylim(-0.5, 0.4)

fig.show()
import cubehelix
cmapper = lambda: cubehelix.cmap(startHue=240,endHue=-300,minSat=1,maxSat=2.5,minLight=.3,maxLight=.8,gamma=.9)

from starfisher.sfhplot import ChiTriptykPlot
fig = plt.figure(figsize=(7, 5))
ctp = ChiTriptykPlot(sfh.full_chi_path, 1, ir_cmd.x_span, ir_cmd.y_span,
                     ir_cmd.dpix, ir_cmd.x_label, ir_cmd.y_label,
                    flipy=True)
ax_obs, ax_mod, ax_chi = ctp.setup_axes(fig)
ctp.plot_obs_in_ax(ax_obs, cmap=cmapper())
ctp.plot_mod_in_ax(ax_mod, cmap=cmapper())
ctp.plot_chi_in_ax(ax_chi, cmap=cubehelix.cmap())
ax_obs.text(0.0, 1.01, "Observed", transform=ax_obs.transAxes, size=8, ha='left')
ax_mod.text(0.0, 1.01, "Model", transform=ax_mod.transAxes, size=8, ha='left')
ax_chi.text(0.0, 1.01, r"$\log \chi^2$", transform=ax_chi.transAxes, size=8, ha='left')
fig.show()<jupyter_output>{'origin': 'lower', 'cmap': <matplotlib.colors.LinearSegmentedColormap object at 0x11752a6d0>, 'extent': [-1.0, 2.0, 25.0, 16.0], 'aspect': 'auto', 'interpolation': 'none'}
{'origin': 'lower', 'cmap': <matplotlib.colors.LinearSegmentedColormap object at 0x116641e90>, 'extent': [-1.0, 2.0, 25.0, 16.0], 'aspect': 'auto', 'interpolation': 'none'}
{'origin': 'lower', 'cmap': <matplotlib.colors.LinearSegmentedColormap object at 0x10e51f2d0>, 'extent': [-1.0, 2.0, 25.0, 16.0], 'aspect': 'auto', 'interpolation': 'none'}
{'fraction_non_alphanumeric': 0.1346704871060172, 'fraction_numerical': 0.06148519579751671, 'mean_word_length': 3.572598253275109, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '260432', 'license': 'MIT', 'n_tokens_mistral': 3587, 'n_tokens_neox': 3229, 'n_words': 593, 'path': 'notebooks/Brick 23 IR.ipynb', 'repo_name': 'jonathansick/androcmd', 'size': '8376'}