Skip to content

Commit

Permalink
Add proto prob
Browse files Browse the repository at this point in the history
  • Loading branch information
@autocorr
autocorr committed Feb 25, 2015
1 parent fb2c6fb commit d364095
Show file tree
Hide file tree
Showing 2 changed files with 90 additions and 31 deletions.
41 changes: 30 additions & 11 deletions besl/bplot/sf_frac.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,21 @@ def get_data():
class ObsData(object):
winr = 25

def __init__(self, col, xlabel):
def __init__(self, col, xlabel, bool_proto=False):
# labels
self.col = col
self.xlabel = xlabel
self.sf_fcol = 'sf_frac_' + col
self.sf_tcol = 'sf_tot_' + col
self.sf_wcol = 'sf_win_' + col
# use proto_prob or boolean probability
if bool_proto:
self.proto_col = 'is_proto'
else:
self.proto_col = 'proto_prob'
# data
self.df = self.get_df()
self.mean_frac = self.df.is_proto.mean()
self.cumix = self.df.index
self.bins = self.df[col].values
self.nbins = self.bins.shape[0]
Expand All @@ -50,7 +56,7 @@ def __init__(self, col, xlabel):
def get_df(self):
df = get_data()
df = df.query('10 < glon_peak < 65')
df = df[df[self.col].notnull()].sort(self.col)
df = df[df[self.col].notnull() & (df[self.col] > 0)].sort(self.col)
stages, labels = dpdf_calc.evo_stages(bgps=df)
df['is_proto'] = False
df.loc[stages[1].index, 'is_proto'] = True
Expand All @@ -60,7 +66,7 @@ def get_df(self):

def cumfrac(self, cumix):
for nn, ii in enumerate(cumix):
nproto = self.df.loc[:ii, 'is_proto'].sum()
nproto = self.df.loc[:ii, self.proto_col].sum()
self.df.loc[ii, self.sf_fcol] = nproto / nn
self.df.loc[ii, self.sf_tcol] = nproto / self.nbins
self.df.loc[ii, self.sf_wcol] = self.window(nn, cumix)
Expand All @@ -73,10 +79,10 @@ def window(self, nn, cumix):
if nn < self.winr:
win = cumix.values[:self.winr+1]
elif nn > self.nbins - self.winr:
win = cumix.values[nn-self.winr:]
win = cumix.values[-self.winr:]
else:
win = cumix.values[nn-self.winr:nn+self.winr+1]
return self.df.loc[win, 'is_proto'].mean()
return self.df.loc[win, self.proto_col].mean()


class ObsPlot(object):
Expand All @@ -88,7 +94,7 @@ def make_fig(self):
ax.set_xscale('log')
ax.set_xlabel(self.od.xlabel)
ax.set_xlim(self.od.xmin, self.od.xmax)
ax.set_ylim(0, 1.05)
ax.set_ylim(0, 1.1)
plt.subplots_adjust(bottom=0.22)
return fig, ax

Expand All @@ -110,20 +116,25 @@ def plot_tot(self):

def plot_win(self):
fig, ax = self.make_fig()
ax.set_ylabel(r'$R_{\rm proto}$' + r' \ ${\rm Win}=' +
str(2 * self.od.winr) + '$' )
ax.set_ylabel(r'$R_{\rm proto}$')
ax.annotate(r'${\rm Window} = ' + str(2 * self.od.winr) + '$',
xy=(0.725, 0.1), xycoords='axes fraction', fontsize=13)
ax.hlines(1, self.od.xmin, self.od.xmax, linestyles='dashed',
colors='0.5')
ax.hlines(self.od.mean_frac, self.od.xmin, self.od.xmax,
linestyles='dotted', colors='0.5')
ax.vlines(self.od.bins, 1.02, 1.07, colors='black', linestyles='solid',
linewidth=0.1)
ax.plot(self.od.bins, self.od.wvals, 'k-', drawstyle='steps')
ax.plot(self.od.bins[:self.od.winr], self.od.wvals[:self.od.winr],
'r-', drawstyle='steps')
color='0.5', linestyle='solid', drawstyle='steps')
ax.plot(self.od.bins[-self.od.winr:], self.od.wvals[-self.od.winr:],
'r-', drawstyle='steps')
color='0.5', linestyle='solid', drawstyle='steps')
util.savefig('{0}_{1}'.format(self.od.sf_wcol, self.od.winr),
close=True)


def plot_window_sizes(wmin=5, wmax=60, step=1):
def plot_window_sizes(wmin=10, wmax=25, step=5):
for ii in range(wmin, wmax+1, step):
print ii
ObsData.winr = ii
Expand Down Expand Up @@ -170,6 +181,14 @@ def plot_all_obs():
op.plot_win()


class MargData(object):
pass


class MargPlot(object):
pass


#####
# Old code to calculate fractions based on molecular detection.
#####
Expand Down
80 changes: 60 additions & 20 deletions besl/clump_match.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@
import numpy as _np
import pandas as _pd
import catalog
from .catalog import (read_bgps, read_cat)
from .image import sample_bgps_img


Expand Down Expand Up @@ -567,7 +566,7 @@ def __init__(self, data):
self.data = data
# BGPS data
self.cat['in_bgps'] = _np.nan
self.bgps = read_bgps(v=self.v).set_index(self.cnum_col)
self.bgps = catalog.read_bgps(v=self.v).set_index(self.cnum_col)
# Process and match
self._add_new_cols()
self._make_haystack()
Expand Down Expand Up @@ -703,7 +702,7 @@ def __init__(self):

def _merge(self):
print '-- Merging data'
merged_data = read_cat('bgps_v210').set_index('v210cnum')
merged_data = catalog.read_cat('bgps_v210').set_index('v210cnum')
for data in self.all_data:
merge_cat = data.matcher.bgps_culled
merged_data = merged_data.merge(merge_cat,
Expand Down Expand Up @@ -796,7 +795,7 @@ class WaterGbt(Data):
def __init__(self):
# Catalog parameters
self.name = 'h2o_gbt'
self.cat = read_cat('gbt_h2o')
self.cat = catalog.read_cat('gbt_h2o')
self.lon_col = 'h2o_glon'
self.lat_col = 'h2o_glat'
self.det_col = 'h2o_f'
Expand All @@ -809,7 +808,7 @@ class WaterArcetri(Data):
def __init__(self):
# Catalog parameters
self.name = 'h2o_arc'
self.cat = read_cat('valdettaro01_arcetri')
self.cat = catalog.read_cat('valdettaro01_arcetri')
self.lon_col = '_Glon'
self.lat_col = '_Glat'
self.det_col = 'h2o_f'
Expand All @@ -822,7 +821,7 @@ class WaterHops(Data):
def __init__(self):
# Catalog parameters
self.name = 'h2o_hops'
self.cat = read_cat('walsh11_hops_h2o')
self.cat = catalog.read_cat('walsh11_hops_h2o')
self.lon_col = 'lPeak_deg'
self.lat_col = 'bPeak_deg'
self.det_col = None
Expand All @@ -835,7 +834,7 @@ class WaterRms(Data):
def __init__(self):
# Catalog parameters
self.name = 'h2o_rms'
self.cat = read_cat('urquhart11_red_msx_h2o')
self.cat = catalog.read_cat('urquhart11_red_msx_h2o')
self.lon_col = '_Glon_1_'
self.lat_col = '_Glat_1_'
self.det_col = 'H2O_1_'
Expand All @@ -848,7 +847,7 @@ class Cornish(Data):
def __init__(self):
# Catalog parameters
self.name = 'corn'
self.cat = read_cat('cornish_uchii')
self.cat = catalog.read_cat('cornish_uchii')
self.lon_col = 'l_deg'
self.lat_col = 'b_deg'
self.det_col = None
Expand All @@ -861,7 +860,7 @@ class Egos(Data):
def __init__(self):
# Catalog parameters
self.name = 'ego'
self.cat = read_cat('ego_all')
self.cat = catalog.read_cat('ego_all')
self.lon_col = '_Glon'
self.lat_col = '_Glat'
self.det_col = None
Expand All @@ -874,7 +873,7 @@ class AmmoniaGbt(Data):
def __init__(self):
# Catalog parameters
self.name = 'nh3_gbt'
self.cat = read_cat('gbt_nh3')
self.cat = catalog.read_cat('gbt_nh3')
self.lon_col = 'glon'
self.lat_col = 'glat'
self.det_col = None
Expand All @@ -887,7 +886,7 @@ class MethoPandian(Data):
def __init__(self):
# Catalog parameters
self.name = 'ch3oh_pandi'
self.cat = read_cat('pandian11')
self.cat = catalog.read_cat('pandian11')
self.lon_col = 'glon'
self.lat_col = 'glat'
self.det_col = None
Expand All @@ -900,7 +899,7 @@ class MethoPestalozzi(Data):
def __init__(self):
# Catalog parameters
self.name = 'ch3oh_pesta'
self.cat = read_cat('pestalozzi05')
self.cat = catalog.read_cat('pestalozzi05')
self.lon_col = '_Glon'
self.lat_col = '_Glat'
self.det_col = None
Expand All @@ -913,7 +912,7 @@ class MethoMmb(Data):
def __init__(self):
# Catalog parameters
self.name = 'ch3oh_mmb'
self.cat = read_cat('mmb_all')
self.cat = catalog.read_cat('mmb_all')
self.lon_col = 'glon'
self.lat_col = 'glat'
self.det_col = None
Expand All @@ -926,7 +925,7 @@ class Higal70(Data):
def __init__(self):
# Catalog parameters
self.name = 'higal70'
self.cat = read_cat('higal_70_clean')
self.cat = catalog.read_cat('higal_70_clean')
self.lon_col = 'GLON'
self.lat_col = 'GLAT'
self.det_col = None
Expand All @@ -939,7 +938,7 @@ class RedSpitzer(Data):
def __init__(self):
# Catalog parameters
self.name = 'robit'
self.cat = read_cat('robitaille08_red_spitzer')
self.cat = catalog.read_cat('robitaille08_red_spitzer')
self.lon_col = '_Glon'
self.lat_col = '_Glat'
self.det_col = 'Class'
Expand All @@ -952,7 +951,7 @@ class RedMsx(Data):
def __init__(self):
# Catalog parameters
self.name = 'red_msx'
self.cat = read_cat('lumsden13_red_msx')
self.cat = catalog.read_cat('lumsden13_red_msx')
self.lon_col = 'glon'
self.lat_col = 'glat'
self.det_col = 'Type'
Expand All @@ -965,7 +964,7 @@ class Molcat(Data):
def __init__(self):
# Catalog parameters
self.name = 'mol'
self.cat = read_cat('shirley13_molcat')
self.cat = catalog.read_cat('shirley13_molcat')
self.lon_col = 'hht_glon'
self.lat_col = 'hht_glat'
self.det_col = 'mol_vlsr_f'
Expand All @@ -978,7 +977,7 @@ class WienenNh3(Data):
def __init__(self):
# Catalog parameters
self.name = 'wien'
cat = read_cat('wienen12_nh3')
cat = catalog.read_cat('wienen12_nh3')
cat = cat[cat['tkin'].notnull()]
self.cat = cat
self.lon_col = '_Glon'
Expand All @@ -993,7 +992,7 @@ class MipsgalCatalog(Data):
def __init__(self):
# Catalog parameters
self.name = 'mipsc'
self.cat = read_cat('mipsgal_catalog_lclip')
self.cat = catalog.read_cat('mipsgal_catalog_lclip')
self.lon_col = 'l'
self.lat_col = 'b'
self.det_col = None
Expand All @@ -1006,7 +1005,7 @@ class MipsgalArchive(Data):
def __init__(self):
# Catalog parameters
self.name = 'mipsa'
self.cat = read_cat('mipsgal_archive_lclip')
self.cat = catalog.read_cat('mipsgal_archive_lclip')
self.lon_col = 'l'
self.lat_col = 'b'
self.det_col = None
Expand All @@ -1015,6 +1014,47 @@ def __init__(self):
self.noise_col = None


###############################################################################
# Protostellar Probability
###############################################################################


class ProtoProb(object):
pyso = 0.5
pagb = 0.4
pcol = 'proto_prob'

def __init__(self, df=None):
if df is None:
self.df = catalog.read_cat('bgps_v210_evo').set_index('v210cnum')
else:
self.df = df
self.df[self.pcol] = 0.0
self.calc()

def calc(self):
# Compute probabilities for R08
rix = self.df.query('robit_n > 0').index
for ii in rix:
nagb = self.df.loc[ii, 'robit_n'] - self.df.loc[ii, 'robit_f']
nyso = self.df.loc[ii, 'robit_f']
self.df.loc[ii, self.pcol] = self.proto_prob(nagb, nyso)
# Correct for certain sources
oix = self.df.query('hg70_eye_f in [1,2,4,5] |'
'ego_n > 0 | '
'red_msx_f > 0 | '
'corn_n > 0 | '
'h2o_f > 0 | '
'ch3oh_f > 0').index
self.df.loc[oix, self.pcol] = 1.0

def proto_prob(self, nagb, nyso):
# Binomial for compliment of zero successes
bagb = 1 - (1 - self.pagb)**nagb
byso = 1 - (1 - self.pyso)**nyso
return byso + bagb - byso * bagb


###############################################################################
# Randomized Cross Matching
###############################################################################
Expand Down

0 comments on commit d364095

Please sign in to comment.