``````

In [1]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import fmin_cg, minimize
import h5py
c = 2.99792458e8   # m/s

``````
``````

In [2]:

def doppler(v):
frac = (1. - v/c) / (1. + v/c)
return np.sqrt(frac)

def gamma(v):
return 1. / np.sqrt(1. - (v/c) ** 2)

def dlndopplerdv(v):
dv = doppler(v)
return -1. * gamma(v) / (2. * c) * (1. / dv  + dv)

def state(v, xs, xps):
'''
outputs: (M, Mp, v, xs, ms, mps, ehs, bees, seas)
M and Mp are the lengths of data and model wavelength grids
v is the RV
xs is the wavelength values of the data grid
ms is the data index m at which there is an interpolated model value
mps is the model index m' from which we interpolate to ms
ehs, bees, and seas go into the coefficients of interpolation
'''
# every input must be 1-d
M = len(xs)
Mp = len(xps)
xps_shifted = xps + np.log(doppler(v))
ms = np.arange(M)
mps = np.searchsorted(xps_shifted, xs, side='left')
good = (mps > 0) * (mps < Mp)
ms = ms[good]
mps = mps[good]
ehs = xps_shifted[mps] - xs[ms]
bees = xs[ms] - xps_shifted[mps - 1]
seas = ehs + bees
return (M, Mp, v, xs, ms, mps, ehs, bees, seas)

def Pdot(state, vec):
# takes state and model flux vector, returns (shifted) model interpolated into data space
# unpack state
M, Mp, v, xs, ms, mps, ehs, bees, seas = state
# do shit
result = np.zeros(M)
result[ms] = vec[mps - 1] * ehs / seas + vec[mps] * bees / seas
return result

def dotP(state, vec):
# takes state and data flux vector, returns data interpolated into (shifted) model space
# unpack state
M, Mp, v, xs, ms, mps, ehs, bees, seas = state
# do shit
result = np.zeros(Mp)
result[mps - 1] += vec[ms] * ehs / seas
result[mps] += vec[ms] * bees / seas
return result

def dotdPdv(state, vec):
# unpack state
M, Mp, v, xs, ms, mps, ehs, bees, seas = state
# do shit
result = np.zeros(Mp)
foos = vec[ms] / seas * dlndopplerdv(v) # * xs[ms] ??
result[mps - 1] += foos
result[mps] -= foos
return result

def dPdotdv(state, vec):
# unpack state
M, Mp, v, xs, ms, mps, ehs, bees, seas = state
# do shit
result = np.zeros(M)
result[ms] = (vec[mps - 1] - vec[mps]) * dlndopplerdv(v) / seas
return result

``````
``````

In [3]:

f = h5py.File('../data/hip54287_15A.hdf5', 'r')

N = 75
data = np.copy(f['data'])[:N,:]
data_xs = np.log(np.copy(f['xs']))
ivars = np.copy(f['ivars'])[:N,:]
true_rvs = np.copy(f['true_rvs'])[:N]
bervs = np.copy(f['berv'])[:N] * -1.e3

for i in xrange(len(data)):
data[i] /= np.median(data[i])

data = np.log(data)

``````
``````

In [4]:

def make_template(all_data, rvs, xs, dx):
"""
`all_data`: `[N, M]` array of pixels
`rvs`: `[N]` array of RVs
`xs`: `[M]` array of wavelength values
`dx`: linear spacing desired for template wavelength grid (A)
"""
(N,M) = np.shape(all_data)
all_xs = np.empty_like(all_data)
for i in range(N):
all_xs[i,:] = xs - np.log(doppler(rvs[i])) # shift to rest frame
all_data, all_xs = np.ravel(all_data), np.ravel(all_xs)
tiny = 10.
template_xs = np.arange(min(all_xs)-tiny*dx, max(all_xs)+tiny*dx, dx)
template_ys = np.nan + np.zeros_like(template_xs)
for i,t in enumerate(template_xs):
ind = (all_xs >= t-dx/2.) & (all_xs < t+dx/2.)
if np.sum(ind) > 0:
template_ys[i] = np.nanmedian(all_data[ind])
ind_nan = np.isnan(template_ys)
template_ys.flat[ind_nan] = np.interp(template_xs[ind_nan], template_xs[~ind_nan], template_ys[~ind_nan]) #np.interp(template_xs[ind_nan], template_xs[~ind_nan], template_ys[~ind_nan])
return template_xs, template_ys

def subtract_template(data_xs, data, model_xs_t, model_ys_t, rvs_t):
(N,M) = np.shape(data)
data_sub = np.copy(data)
for n,v in enumerate(rvs_t):
s = state(v, data_xs, model_xs_t)
model_ys_t_shifted = Pdot(s, model_ys_t)
data_sub[n,:] -= np.ravel(model_ys_t_shifted)
if n == 0:
plt.plot(data_xs, data[n,:], color='k')
plt.plot(data_xs, data_sub[n,:], color='blue')
plt.plot(data_xs, np.ravel(model_ys_t_shifted), color='red')
return data_sub

``````
``````

In [5]:

x0_star = -np.copy(bervs)
x0_star -= np.mean(x0_star)
x0_t = np.zeros(N)
model_xs_star, model_ys_star = make_template(data, x0_star, data_xs, np.log(6000.01) - np.log(6000.))
model_xs_t, model_ys_t = make_template(data, x0_t, data_xs, np.log(6000.01) - np.log(6000.))

``````
``````

In [6]:

plt.plot(model_xs_star, model_ys_star, color='k')
plt.plot(model_xs_t, model_ys_t, color='red')
#plt.plot(data_xs, data[0,:], color='blue')

``````
``````

Out[6]:

[<matplotlib.lines.Line2D at 0x112def1d0>]

``````
``````

In [19]:

plt.hist(model_ys_star, color='k', alpha=0.5, bins=1000)
plt.hist(model_ys_t, color='red', alpha=0.5, bins=1000)
plt.yscale('log', nonposy='clip')
#plt.ylim((0,100))

``````
``````

``````
``````

In [20]:

stepsize = 0.001
for n in range(75):
tmp = np.copy(x0_star)
s = state(tmp[n], data_xs, model_xs_star)
ddotP1 = dotdPdv(s, data[n,:])
tmp[n] += stepsize
s = state(tmp[n], data_xs, model_xs_star)
dotP0 = dotP(s, data[n,:])
tmp[n] -= 2. * stepsize
s = state(tmp[n], data_xs, model_xs_star)
dotP2 = dotP(s, data[n,:])
print ddotP1[1000:1004]
print (dotP0 - dotP2)[1000:1004]/(2 * stepsize)

``````
``````

[ -2.69414869e-05  -6.94456849e-06   2.71558932e-05   1.91587986e-05]
[ -2.69371981e-05  -6.94346298e-06   2.71515703e-05   1.91557487e-05]
[  1.08835211e-05  -1.58762189e-05  -3.85632344e-05   1.49649364e-05]
[  1.08846866e-05  -1.58779190e-05  -3.85673639e-05   1.49665388e-05]
[ -1.09399898e-05  -1.39728627e-05  -7.35531635e-06  -2.70795854e-06]
[ -1.09382482e-05  -1.39706383e-05  -7.35414546e-06  -2.70752746e-06]
[  1.21707276e-05  -1.82132346e-05   2.36035890e-06   1.03890811e-05]
[  1.21720309e-05  -1.82151849e-05   2.36061165e-06   1.03901935e-05]
[  4.69326451e-06  -7.87548860e-07  -5.49822535e-06  -5.91914654e-06]
[  4.69251741e-06  -7.87423494e-07  -5.49735011e-06  -5.91820430e-06]
[ -1.11401625e-05  -2.86261252e-05  -4.50569516e-06   2.34632785e-05]
[ -1.11383891e-05  -2.86215683e-05  -4.50497792e-06   2.34595435e-05]
[ -6.86534588e-06  -1.55065481e-05  -2.10916303e-05  -2.21527018e-05]
[ -6.86608104e-06  -1.55082086e-05  -2.10938889e-05  -2.21550740e-05]
[  1.97399748e-06  -9.10895542e-06  -1.90056844e-05  -4.92416786e-06]
[  1.97420886e-06  -9.10993083e-06  -1.90077196e-05  -4.92469515e-06]
[  1.02686893e-06   8.50856886e-07  -3.30267066e-06  -1.11657338e-05]
[  1.02697889e-06   8.50947997e-07  -3.30302432e-06  -1.11669294e-05]
[  3.63395432e-06   5.21012648e-06  -4.09049304e-06  -6.87671910e-06]
[  3.63434345e-06   5.21068438e-06  -4.09093105e-06  -6.87745546e-06]
[ -7.35780621e-06  -7.26975974e-06  -4.68216470e-06   7.17014069e-07]
[ -7.35859411e-06  -7.27053821e-06  -4.68266608e-06   7.17090849e-07]
[  3.50817919e-06   1.10866356e-05   2.14282691e-06  -9.79585836e-06]
[  3.50762074e-06   1.10848708e-05   2.14248581e-06  -9.79429899e-06]
[ -7.78248731e-06   2.38362274e-06   7.66367864e-06  -3.95597969e-06]
[ -7.78332068e-06   2.38387799e-06   7.66449929e-06  -3.95640330e-06]
[ -2.14377645e-05  -2.99661646e-06   1.57272591e-05   2.72940309e-06]
[ -2.14400601e-05  -2.99693734e-06   1.57289432e-05   2.72969536e-06]
[ -6.38497369e-06   1.12720687e-05  -6.01597850e-06  -1.46998457e-05]
[ -6.38565742e-06   1.12732757e-05  -6.01662271e-06  -1.47014198e-05]
[  1.51142504e-05   1.68079777e-07  -3.62598621e-06   2.66228613e-06]
[  1.51118444e-05   1.68053021e-07  -3.62540900e-06   2.66186233e-06]
[  2.87591355e-06   9.42579207e-06   1.83614937e-05   1.68586146e-05]
[  2.87622151e-06   9.42680141e-06   1.83634599e-05   1.68604199e-05]
[ -9.69233040e-06  -1.88281730e-05  -7.36464953e-06   1.65870040e-06]
[ -9.69336828e-06  -1.88301892e-05  -7.36543816e-06   1.65887802e-06]
[ -7.93818103e-06   7.60044761e-07   2.36469402e-06   7.84221832e-07]
[ -7.93903107e-06   7.60126149e-07   2.36494723e-06   7.84305808e-07]
[ -1.87995472e-06  -5.87892888e-06   3.05226368e-06  -2.88416463e-06]
[ -1.88015603e-06  -5.87955841e-06   3.05259052e-06  -2.88447347e-06]
[ -6.68360361e-07   2.61492264e-06   1.63649370e-06  -4.71548303e-06]
[ -6.68253966e-07   2.61450638e-06   1.63623319e-06  -4.71473239e-06]
[  2.17332166e-05  -4.40524014e-06  -1.87323165e-05   2.18437740e-05]
[  2.17355437e-05  -4.40571185e-06  -1.87343224e-05   2.18461130e-05]
[  1.12603414e-07  -1.61443041e-05  -3.98039658e-06   1.16053215e-05]
[  1.12615472e-07  -1.61460329e-05  -3.98082281e-06   1.16065642e-05]
[ -4.66331285e-06   7.00422259e-07   1.92289612e-06   6.13962564e-06]
[ -4.66381221e-06   7.00497262e-07   1.92310203e-06   6.14028309e-06]
[ -2.45920259e-06   1.04585268e-06   4.14684548e-07  -8.08942342e-06]
[ -2.45881111e-06   1.04568619e-06   4.14618534e-07  -8.08813569e-06]
[ -1.28553361e-05   3.28373827e-05   1.48459233e-05  -1.07843380e-05]
[ -1.28567127e-05   3.28408990e-05   1.48475130e-05  -1.07854928e-05]
[ -1.43629653e-05  -9.48200965e-06  -1.50581636e-05  -1.04056967e-05]
[ -1.43645033e-05  -9.48302500e-06  -1.50597761e-05  -1.04068110e-05]
[  3.61584822e-06  -8.93289922e-06  -1.78327179e-05  -4.31986828e-06]
[  3.61623542e-06  -8.93385578e-06  -1.78346275e-05  -4.32033086e-06]
[ -1.08515091e-05  -1.09732874e-05  -1.51991889e-06   5.82689063e-06]
[ -1.08497817e-05  -1.09715406e-05  -1.51967693e-06   5.82596306e-06]
[ -8.41515188e-06  -3.33616029e-06   3.78660580e-06  -3.53795613e-06]
[ -8.41381228e-06  -3.33562921e-06   3.78600301e-06  -3.53739293e-06]
[ -1.23480262e-05  -1.01245185e-05  -4.44659982e-06   6.48380210e-06]
[ -1.23493484e-05  -1.01256027e-05  -4.44707597e-06   6.48449638e-06]
[  4.84238797e-06  -1.66661862e-06  -1.01406253e-05  -7.41806731e-06]
[  4.84290649e-06  -1.66679708e-06  -1.01417112e-05  -7.41886163e-06]
[  9.58135800e-07  -1.07593735e-05   1.43905885e-06   1.03435788e-05]
[  9.57983276e-07  -1.07576607e-05   1.43882977e-06   1.03419322e-05]
[ -7.09458650e-06  -1.82460501e-05   5.19361040e-06   2.23916172e-05]
[ -7.09345712e-06  -1.82431456e-05   5.19278363e-06   2.23880527e-05]
[  7.35231739e-06  -9.93467066e-06  -7.98238251e-06   7.16226883e-06]
[  7.35114698e-06  -9.93308916e-06  -7.98111179e-06   7.16112867e-06]
[ -3.27632359e-06  -5.41647441e-06  -2.07983449e-07   5.72209823e-06]
[ -3.27667442e-06  -5.41705441e-06  -2.08005720e-07   5.72271094e-06]
[ -3.11984380e-06  -8.28075044e-06  -4.74368300e-06   7.49334194e-06]
[ -3.12017788e-06  -8.28163716e-06  -4.74419096e-06   7.49414435e-06]
[ -9.94136364e-06   7.45552291e-06   4.84384865e-06  -9.45410368e-06]
[ -9.94242819e-06   7.45632127e-06   4.84436734e-06  -9.45511605e-06]
[  1.62835534e-06  -1.36186041e-05  -4.30741735e-06   4.39338963e-06]
[  1.62809613e-06  -1.36164362e-05  -4.30673167e-06   4.39269026e-06]
[  6.41251502e-06  -6.13140043e-06  -1.37005751e-05  -1.38118885e-05]
[  6.41320170e-06  -6.13205700e-06  -1.37020422e-05  -1.38133675e-05]
[  1.70347672e-05  -2.50173899e-06  -1.50537772e-05  -1.84693475e-05]
[  1.70320555e-05  -2.50134074e-06  -1.50513808e-05  -1.84664074e-05]
[ -7.04603948e-06  -4.86002943e-06   2.19087974e-06   1.00467452e-05]
[ -7.04679400e-06  -4.86054986e-06   2.19111434e-06   1.00478210e-05]
[ -1.36706890e-05   2.64167499e-06   1.57797743e-06  -5.08750344e-06]
[ -1.36721529e-05   2.64195787e-06   1.57814641e-06  -5.08804823e-06]
[  1.39555429e-05  -2.12941131e-05  -6.27970099e-06   9.67860545e-06]
[  1.39570373e-05  -2.12963934e-05  -6.28037344e-06   9.67964187e-06]
[ -1.19524170e-05  -1.54647905e-05  -8.24720832e-06   1.46216942e-05]
[ -1.19536969e-05  -1.54664465e-05  -8.24809146e-06   1.46232600e-05]
[ -4.27678647e-06  -9.85969028e-06  -1.26821082e-05  -9.56557463e-06]
[ -4.27724444e-06  -9.86074609e-06  -1.26834662e-05  -9.56659895e-06]
[ -1.40098128e-05   7.26369953e-06  -3.16579752e-06   6.25296278e-06]
[ -1.40113130e-05   7.26447735e-06  -3.16613653e-06   6.25363237e-06]
[ -9.97560743e-06  -4.09591161e-06   8.22595960e-06  -9.68507827e-07]
[ -9.97401940e-06  -4.09525958e-06   8.22465010e-06  -9.68353649e-07]
[ -8.71614907e-06   4.96409546e-06   1.03882100e-05  -1.16243092e-05]
[ -8.71476154e-06   4.96330523e-06   1.03865563e-05  -1.16224587e-05]
[ -2.35271812e-05  -9.29676794e-06   1.08432397e-05   1.40919279e-05]
[ -2.35297004e-05  -9.29776343e-06   1.08444008e-05   1.40934368e-05]
[ -1.11442207e-05  -5.74944414e-06   2.68308843e-06  -1.26846300e-06]
[ -1.11424467e-05  -5.74852888e-06   2.68266130e-06  -1.26826107e-06]
[ -1.07824925e-05   4.53453999e-07   1.19884002e-05   1.58461086e-06]
[ -1.07836470e-05   4.53502553e-07   1.19896839e-05   1.58478053e-06]
[ -1.25538086e-05  -4.23916055e-07   1.26340656e-05  -4.03983433e-06]
[ -1.25551529e-05  -4.23961448e-07   1.26354185e-05  -4.04026691e-06]
[  3.42032851e-06  -2.92837622e-06  -8.23794195e-06  -3.02597077e-06]
[  3.41978403e-06  -2.92791005e-06  -8.23663054e-06  -3.02548906e-06]
[ -1.11338970e-05  -1.61409619e-06   7.30626143e-06  -4.45782986e-07]
[ -1.11321246e-05  -1.61383924e-06   7.30509835e-06  -4.45712022e-07]
[ -1.45816229e-05   3.23552170e-07   1.63568820e-05  -2.99082923e-06]
[ -1.45793017e-05   3.23500664e-07   1.63542781e-05  -2.99035312e-06]
[  7.30132061e-06   1.80006838e-06  -9.84225850e-06  -3.77109393e-06]
[  7.30210243e-06   1.80026113e-06  -9.84331241e-06  -3.77149773e-06]
[ -5.12447870e-07   5.99261569e-07  -5.67783860e-06  -7.63889891e-06]
[ -5.12502742e-07   5.99325737e-07  -5.67844658e-06  -7.63971688e-06]
[ -1.06938104e-05  -1.71531663e-05  -3.67171650e-06   5.05683255e-07]
[ -1.06949555e-05  -1.71550031e-05  -3.67210967e-06   5.05737404e-07]
[ -4.62454525e-06  -1.07726123e-05  -5.61091048e-06  -4.53259281e-06]
[ -4.62380908e-06  -1.07708974e-05  -5.61001730e-06  -4.53187128e-06]
[ -7.75065215e-06   5.35809536e-06   9.18357452e-06   3.30155722e-06]
[ -7.74941835e-06   5.35724242e-06   9.18211262e-06   3.30103166e-06]
[ -1.21339547e-05  -6.04973896e-06   2.34628188e-05  -1.96544161e-05]
[ -1.21352541e-05  -6.05038677e-06   2.34653312e-05  -1.96565208e-05]
[ -7.37879939e-06  -5.10784891e-06  -3.44242803e-06  -3.35431757e-06]
[ -7.37958950e-06  -5.10839585e-06  -3.44279664e-06  -3.35467675e-06]
[ -6.01279611e-06  -2.18412037e-06  -2.64615917e-07  -1.52454619e-06]
[ -6.01183893e-06  -2.18377267e-06  -2.64573793e-07  -1.52430350e-06]
[ -2.65926139e-06   2.30690110e-07   5.31997677e-07  -2.80109006e-06]
[ -2.65954614e-06   2.30714813e-07   5.32054643e-07  -2.80139000e-06]
[  2.06359474e-06   2.86823513e-06  -5.30392454e-06  -5.78120749e-06]
[  2.06381570e-06   2.86854226e-06  -5.30449248e-06  -5.78182654e-06]
[  7.22113353e-06   2.51857843e-06  -1.44356877e-05  -3.65012506e-06]
[  7.22190676e-06   2.51884812e-06  -1.44372334e-05  -3.65051591e-06]
[ -2.03366766e-05  -1.19876490e-05   6.66761342e-06   8.00809099e-06]
[ -2.03388543e-05  -1.19889327e-05   6.66832738e-06   8.00894849e-06]
[  3.61452018e-06   8.03406129e-06   5.38371132e-06  -1.60054390e-06]
[  3.61490722e-06   8.03492156e-06   5.38428781e-06  -1.60071528e-06]
[ -1.49739619e-05  -4.00309627e-07   1.76928784e-05   2.41608535e-06]
[ -1.49715782e-05  -4.00245902e-07   1.76900619e-05   2.41570073e-06]
[ -2.42447227e-05  -1.75801984e-05   1.76034960e-05   7.98550956e-06]
[ -2.42473189e-05  -1.75820809e-05   1.76053809e-05   7.98636464e-06]
[  1.95747883e-06  -3.49703430e-06  -9.15688908e-06   3.49124224e-07]
[  1.95768844e-06  -3.49740876e-06  -9.15786959e-06   3.49161606e-07]
[ -3.44361158e-05  -6.97518445e-06   9.31681853e-06   4.15972807e-06]
[ -3.44306340e-05  -6.97407409e-06   9.31533541e-06   4.15906589e-06]
[ -1.24975175e-05  -1.13488162e-05  -6.48332609e-06  -3.11956780e-06]
[ -1.24988557e-05  -1.13500314e-05  -6.48402034e-06  -3.11990185e-06]
[ -7.48686927e-06  -1.11808444e-05  -6.57873073e-06  -6.63933637e-06]
[ -7.48567746e-06  -1.11790645e-05  -6.57768348e-06  -6.63827948e-06]

``````
``````

In [21]:

stepsize = 0.01
rv_steps = np.arange(0, 10, stepsize)
M = len(data[0,:])
Pdots = np.zeros((len(rv_steps), M))
dPdots = np.zeros((len(rv_steps), M))
for i in xrange(len(rv_steps)):
drv = rv_steps[i]
s = state(drv, data_xs, model_xs_star)
Pdots[i,:] = Pdot(s, model_ys_star)
dPdots[i,:] = dPdotdv(s, model_ys_star)

print (Pdots[-1] - Pdots[0])[1000:1004]
print (np.sum(dPdots*stepsize, axis=0))[1000:1004]

``````
``````

[  1.57089485e-05  -2.40025100e-05   3.77128655e-05   1.44735039e-05]
[  1.57246733e-05  -2.40265368e-05   3.77506165e-05   1.44879921e-05]

``````
``````

In [22]:

def rv_lnprior(rvs):
return -0.5 * np.mean(rvs)**2/1.**2

def drv_lnprior_dv(rvs):
return np.zeros_like(rvs) - np.mean(rvs)/1.**2/len(rvs)

def lnlike_star(rvs_star, rvs_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t):
try:
N = len(rvs_star)
except:
N = 1
lnlike = 0.
dlnlike_dv = np.zeros(N)
for n in range(N):
state_star = state(rvs_star[n], data_xs, model_xs_star)
pd_star = Pdot(state_star, model_ys_star)
state_t = state(rvs_t[n], data_xs, model_xs_t)
pd_t = Pdot(state_t, model_ys_t)
pd = pd_star + pd_t
lnlike += -0.5 * np.sum((data[n,:] - pd)**2 * ivars[n,:])
dpd_dv = dPdotdv(state_star, model_ys_star)
dlnlike_dv[n] = np.sum((data[n,:] - pd) * ivars[n,:] * dpd_dv)
lnpost = lnlike + rv_lnprior(rvs_star)
dlnpost_dv = dlnlike_dv + drv_lnprior_dv(rvs_star)
return -1 * lnpost, -1 * dlnpost_dv

def lnlike_t(rvs_t, rvs_star, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t):
try:
N = len(rvs_t)
except:
N = 1
lnlike = 0.
dlnlike_dv = np.zeros(N)
for n in range(N):
state_star = state(rvs_star[n], data_xs, model_xs_star)
pd_star = Pdot(state_star, model_ys_star)
state_t = state(rvs_t[n], data_xs, model_xs_t)
pd_t = Pdot(state_t, model_ys_t)
pd = pd_star + pd_t
lnlike += -0.5 * np.sum((data[n,:] - pd)**2 * ivars[n,:])
dpd_dv = dPdotdv(state_t, model_ys_t)
dlnlike_dv[n] = np.sum((data[n,:] - pd) * ivars[n,:] * dpd_dv)
lnpost = lnlike + rv_lnprior(rvs_t)
dlnpost_dv = dlnlike_dv + drv_lnprior_dv(rvs_t)
return -1 * lnpost, -1 * dlnpost_dv

``````
``````

In [23]:

def model_ys_lnprior(w):
return -0.5 * np.sum(w**2)/100.**2

def dmodel_ys_lnprior_dw(w):
return -1.*w / 100.**2

def dlnlike_star_dw_star(model_ys_star, rvs_star, rvs_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t):
try:
N = len(rvs_star)
except:
N = 1
lnlike = 0.
Mp = len(model_xs_star)
dlnlike_dw = np.zeros(Mp)
for n in range(N):
state_star = state(rvs_star[n], data_xs, model_xs_star)
pd_star = Pdot(state_star, model_ys_star)
state_t = state(rvs_t[n], data_xs, model_xs_t)
pd_t = Pdot(state_t, model_ys_t)
pd = pd_star + pd_t
dp_star = dotP(state_star, (data[n,:] - pd)*ivars[n,:])
lnlike += -0.5 * np.sum((data[n,:] - pd)**2 * ivars[n,:])
dlnlike_dw += dp_star
lnprior = model_ys_lnprior(model_ys_star)
dlnprior = dmodel_ys_lnprior_dw(model_ys_star)
return -lnlike - lnprior, -dlnlike_dw - dlnprior

def dlnlike_t_dw_t(model_ys_t, rvs_star, rvs_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t):
try:
N = len(rvs_t)
except:
N = 1
lnlike = 0.
Mp = len(model_xs_t)
dlnlike_dw = np.zeros(Mp)
for n in range(N):
state_star = state(rvs_star[n], data_xs, model_xs_star)
pd_star = Pdot(state_star, model_ys_star)
state_t = state(rvs_t[n], data_xs, model_xs_t)
pd_t = Pdot(state_t, model_ys_t)
pd = pd_star + pd_t
dp_t = dotP(state_t, (data[n,:] - pd)*ivars[n,:])
lnlike += -0.5 * np.sum((data[n,:] - pd)**2 * ivars[n,:])
dlnlike_dw += dp_t
lnprior = model_ys_lnprior(model_ys_t)
dlnprior = dmodel_ys_lnprior_dw(model_ys_t)
return -lnlike - lnprior, -dlnlike_dw - dlnprior

``````
``````

In [24]:

print dlnlike_star_dw_star(model_ys_star, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)

stepsize = 0.001
lnlike1, dlnlike_dv1 = dlnlike_star_dw_star(model_ys_star, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)
for n in range(75):
tmp = np.copy(model_ys_star)
tmp[n] += stepsize
lnlike0, dlnlike_dv0 = dlnlike_star_dw_star(tmp, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)
tmp[n] -= 2. * stepsize
lnlike2, dlnlike_dv2 = dlnlike_star_dw_star(tmp, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)
print dlnlike_dv1[n]
print (lnlike0 - lnlike2)/(2 * stepsize)

``````
``````

(443867.67954574793, array([  9.50767486e-07,   9.50767486e-07,   9.50767486e-07, ...,
6.34706276e-07,   6.34706276e-07,   6.34706276e-07]))
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
9.50767486473e-07
9.60426405072e-07
1223.59356171
1223.59356171
1998.69177424
1998.69177424
2635.63451317
2635.63451316
2205.27875769
2205.27875764
2095.19516647
2095.19516645
1946.15191326
1946.15191329
2236.67956905
2236.679569
3356.07939607
3356.07939609
4339.5678441
4339.56784406
3928.3188656
3928.3188656
3754.36230137
3754.36230141
3834.36230047
3834.36230055
4306.59634869
4306.59634867
4071.24057232
4071.24057235
4276.08536278
4276.08536277
5043.64615078
5043.64615079
4724.93108553
4724.93108551
4522.60381072
4522.60381071
4971.68778455
4971.68778451
6026.41696242
6026.41696241
5923.96609084
5923.96609086
4419.85452147
4419.85452149
4648.7790714
4648.77907143
5861.59717886
5861.59717882
6572.50842804
6572.50842807
6949.82137528
6949.82137528
6907.21436164
6907.21436171
6034.05919743
6034.0591975
5314.3812829
5314.38128286
6113.91988666
6113.91988667
5556.1903321
5556.19033208
5392.96014203
5392.96014202
3637.22492375
3637.22492373
2490.44779373
2490.4477938
3731.77816697
3731.77816704
4864.07182242
4864.07182241
4674.41610953
4674.41610948
6661.48747454
6661.48747454
9138.08251989
9138.08251987
9450.26471714
9450.2647171
9570.32368802
9570.32368801
7275.32227144
7275.32227142
6701.18307392
6701.18307398
7659.67884782
7659.67884782
8529.51249385
8529.51249381
7948.11510088
7948.11510091
6829.42912097
6829.42912093
7544.55827282
7544.55827276
6110.54246736
6110.54246742
5419.93111843
5419.93111841
7456.84305251
7456.84305261
8714.45948216
8714.45948211
9567.60445542
9567.60445543
9863.71110299
9863.71110301
9530.69639671
9530.69639672
9543.58445681
9543.58445684
10130.3686064
10130.3686064
10588.0834648
10588.0834649
10493.4517205
10493.4517205
9969.08964211
9969.08964211
11311.5402837
11311.5402837
12022.9818494
12022.9818493
11508.1566366
11508.1566366
13635.5986136
13635.5986136
12343.9043495
12343.9043495

``````
``````

In [25]:

print dlnlike_t_dw_t(model_ys_t, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t)

stepsize = 0.001
lnlike1, dlnlike_dv1 = dlnlike_t_dw_t(model_ys_t, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t)
for n in range(75):
tmp = np.copy(model_ys_t)
tmp[n] += stepsize
lnlike0, dlnlike_dv0 = dlnlike_t_dw_t(tmp, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t)
tmp[n] -= 2. * stepsize
lnlike2, dlnlike_dv2 = dlnlike_t_dw_t(tmp, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t)
print dlnlike_dv1[n]
print (lnlike0 - lnlike2)/(2 * stepsize)

``````
``````

(443867.67859593226, array([  7.83655022e-07,   7.83655022e-07,   7.83655022e-07, ...,
2.38896489e-07,   2.38896489e-07,   2.38896489e-07]))
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
7.83655022199e-07
7.85803422332e-07
16515.474385
16515.4743849
17126.244434
17126.244434
17937.7633678
17937.7633678
15809.4214444
15809.4214444
15504.5267249
15504.5267249
15873.3700128
15873.3700128
17092.8434363
17092.8434363
19766.0081874
19766.0081873
19708.6072279
19708.6072279
16825.9354598
16825.9354597
16954.5417969
16954.5417968
17645.9807546
17645.9807545
17834.7031069
17834.7031069
14514.6971902
14514.6971902
16652.5931886
16652.5931886
17201.8783134
17201.8783134
15136.3476783
15136.3476783
12545.7949184
12545.7949184
10541.8038299
10541.8038298
11905.0079175
11905.0079176
11474.797834
11474.7978341
10956.9125493
10956.9125492
12652.873218
12652.8732181
14740.4705377
14740.4705377
14535.211032
14535.211032
13442.1561094
13442.1561094
13885.7715629
13885.7715629
14491.9365644
14491.9365644
14645.2945429
14645.2945429
13877.8676951
13877.8676952
12557.0087563
12557.0087562
11862.8349886
11862.8349886
11932.1233193
11932.1233192
13268.0410693
13268.0410693
15012.726007
15012.7260071
15277.0179665
15277.0179665
15023.0962087
15023.0962088
15181.9117483
15181.9117483
14617.5315249
14617.5315251
13449.7657096
13449.7657097
12744.1140771
12744.114077
13362.6712458
13362.6712457
14863.567079
14863.5670791
15214.0988503
15214.0988502
15335.448664
15335.4486641
16625.441894
16625.441894
16929.7418623
16929.7418623
16376.6661863
16376.6661863
16437.7311389
16437.7311389
16706.0355494
16706.0355494
16608.8420722
16608.8420721
16508.3646475
16508.3646474
15370.9107439
15370.9107439
17500.3561716
17500.3561715
18867.6021575
18867.6021576
18480.0834837
18480.0834838
17640.3140992
17640.3140993
17722.3894821
17722.3894821
18181.3407566
18181.3407568
20745.6761186
20745.6761186
21059.9813914
21059.9813915
22273.6111451
22273.611145
21083.8422285
21083.8422285
20911.5311126
20911.5311127
21916.4786535
21916.4786535

``````
``````

In [26]:

dw = dlnlike_star_dw_star(model_ys_star, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)[1]

plt.plot(model_xs_star, dw, 'k')
plt.show()
plt.plot(model_xs_star, model_ys_star, 'k')

``````
``````

Out[26]:

[<matplotlib.lines.Line2D at 0x1167d2750>]

``````
``````

In [27]:

w = np.copy(model_ys_star + np.random.normal(0, 0.01, len(model_ys_star)))
for i in range(50):
lnlike, dlnlike_dw = dlnlike_star_dw_star(w, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)
w -= 5e-7 * dlnlike_dw

``````
``````

570942.838881
366115.363987
353545.678639
351335.272246
350662.484882
350409.027794
350299.525654
350246.005071
350216.746721
350199.180521
350187.840125
350180.114364
350174.638281
350170.638885
350167.648487
350165.369031
350163.602826
350162.214641
350161.109674
350160.22012
350159.496641
350158.902791
350158.411276
350158.001396
350157.657269
350157.366572
350157.119646
350156.908843
350156.728054
350156.572359
350156.437759
350156.320983
350156.219338
350156.130588
350156.05287
350155.984625
350155.924535
350155.871489
350155.82454
350155.782884
350155.745832
350155.71279
350155.683251
350155.656777
350155.632986
350155.611551
350155.592185
350155.574642
350155.558703
350155.544181

``````
``````

In [28]:

plt.plot(model_xs_star, model_ys_star-w, 'k')
plt.show()

plt.plot(model_xs_star, model_ys_star, 'k')
plt.plot(model_xs_star, w, 'r')

``````
``````

Out[28]:

[<matplotlib.lines.Line2D at 0x118b73110>]

``````
``````

In [29]:

def improve_telluric_model(model_ys_t, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t):
w = np.copy(model_ys_t)
for i in range(50):
lnlike, dlnlike_dw = dlnlike_t_dw_t(w, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t)
w -= 5e-7 * dlnlike_dw
return w

def improve_star_model(model_ys_star, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t):
w = np.copy(model_ys_star + np.random.normal(0, 0.01, len(model_ys_star)))
for i in range(50):
lnlike, dlnlike_dw = dlnlike_star_dw_star(w, x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)
w -= 5e-7 * dlnlike_dw
return w

``````
``````

In [30]:

plt.plot(model_xs_t, model_ys_t-w, 'k')
plt.show()

plt.plot(model_xs_t, model_ys_t, 'k')
plt.plot(model_xs_t, w, 'r')

``````
``````

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-30-8fc9aa081c07> in <module>()
----> 1 plt.plot(model_xs_t, model_ys_t-w, 'k')
2 plt.show()
3
4 plt.plot(model_xs_t, model_ys_t, 'k')
5 plt.plot(model_xs_t, w, 'r')

ValueError: operands could not be broadcast together with shapes (1540,) (1625,)

``````
``````

In [ ]:

print lnlike_star(x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)
print lnlike_star(true_rvs, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)

``````
``````

In [ ]:

test_rvs = np.linspace(-1.e3+x0_star[0], 1.e3+x0_star[0], 100)
chisq = np.zeros_like(test_rvs)
for i,v in enumerate(test_rvs):
tmp = np.copy(x0_star)
tmp[0] = v
lnlike0, dlnlike_dv0 = lnlike_star(tmp, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)
chisq[i] = lnlike0

plt.plot(test_rvs, chisq)
plt.axvline(x0_star[0])

``````
``````

In [ ]:

stepsize = 0.001
lnlike1, dlnlike_dv1 = lnlike_star(x0_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)
for n in range(75):
tmp = np.copy(x0_star)
tmp[n] += stepsize
lnlike0, dlnlike_dv0 = lnlike_star(tmp, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)
tmp[n] -= 2. * stepsize
lnlike2, dlnlike_dv2 = lnlike_star(tmp, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)
print dlnlike_dv1[n]
print (lnlike0 - lnlike2)/(2 * stepsize)

``````
``````

In [ ]:

stepsize = 0.001
lnlike1, dlnlike_dv1 = lnlike_t(x0_t, x0_star, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)
for n in range(75):
tmp = np.copy(x0_t)
tmp[n] += stepsize
lnlike0, dlnlike_dv0 = lnlike_t(tmp, x0_star, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)
tmp[n] -= 2. * stepsize
lnlike2, dlnlike_dv2 = lnlike_t(tmp, x0_star, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t)
print dlnlike_dv1[n]
print (lnlike0 - lnlike2)/(2 * stepsize)

``````
``````

In [ ]:

soln_star =  minimize(lnlike_star, x0_star, args=(x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t),
method='BFGS', jac=True, options={'disp':True, 'gtol':1.e-2, 'eps':1.5e-5})['x']
print np.std(soln_star + true_rvs)

``````
``````

In [ ]:

plt.scatter(np.arange(N), x0_star)
plt.scatter(np.arange(N), soln_star)
plt.scatter(np.arange(N), -true_rvs)
plt.show()

``````
``````

In [ ]:

print soln_star - x0_star

``````
``````

In [ ]:

soln_star =  minimize(lnlike_star, x0_star, args=(x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t),
method='BFGS', jac=True, options={'disp':True, 'gtol':1.e-2, 'eps':1.5e-5})['x']
soln_t =  minimize(lnlike_t, x0_t, args=(soln_star, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t),
method='BFGS', jac=True, options={'disp':True, 'gtol':1.e-2, 'eps':1.5e-5})['x']

x0_star = soln_star
x0_t = soln_t
print np.std(x0_star + true_rvs)
print np.std(x0_t)

``````
``````

In [ ]:

plt.scatter(np.arange(N), soln_star + true_rvs)

``````
``````

In [ ]:

data_star = subtract_template(data_xs, data, model_xs_t, model_ys_t, x0_t)

``````
``````

In [ ]:

data_t = subtract_template(data_xs, data, model_xs_star, model_ys_star, x0_star)

``````
``````

In [ ]:

plt.plot(data_xs, data[0,:], color='black')
plt.plot(model_xs_star, model_ys_star, color='red')
plt.plot(model_xs_t, model_ys_t, color='green')

``````
``````

In [ ]:

soln_star =  minimize(lnlike_star, x0_star, args=(x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t),
method='BFGS', jac=True, options={'disp':True, 'gtol':1.e-2, 'eps':1.5e-5})['x']
soln_t =  minimize(lnlike_t, x0_t, args=(soln_star, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t),
method='BFGS', jac=True, options={'disp':True, 'gtol':1.e-2, 'eps':1.5e-5})['x']

print np.std(soln_star + true_rvs)
print np.std(soln_t)

``````
``````

In [ ]:

plt.scatter(np.arange(N), soln_star + true_rvs)

``````
``````

In [ ]:

for n in range(5):
soln_star =  minimize(lnlike_star, x0_star, args=(x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t),
method='BFGS', jac=True, options={'disp':True, 'gtol':1.e-2, 'eps':1.5e-5})['x']

model_ys_t = improve_telluric_model(model_ys_t, soln_star, x0_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t)
model_ys_star = improve_star_model(model_ys_star, soln_star, x0_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)

soln_t =  minimize(lnlike_t, x0_t, args=(soln_star, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t, model_ys_t),
method='BFGS', jac=True, options={'disp':True, 'gtol':1.e-2, 'eps':1.5e-5})['x']

model_ys_t = improve_telluric_model(model_ys_t, soln_star, soln_t, data_xs, data, ivars, model_xs_star, model_ys_star, model_xs_t)
model_ys_star = improve_star_model(model_ys_star, soln_star, soln_t, data_xs, data, ivars, model_xs_star, model_xs_t, model_ys_t)

x0_star = soln_star
x0_t = soln_t

print "iter {0}: star std = {1:.2f}, telluric std = {2:.2f}".format(n, np.std(soln_star + true_rvs), np.std(soln_t))
plt.scatter(np.arange(N), soln_star + true_rvs, color='k')
plt.scatter(np.arange(N), soln_t, color='red')
plt.show()

``````
``````

In [ ]:

plt.scatter(bervs, soln_star+true_rvs)

``````
``````

In [ ]:

plt.plot(data_xs, data[0,:], color='k')
plt.plot(data_xs, data_star[0,:] + data_t[0,:], color='red')

``````
``````

In [ ]:

plt.plot(data_xs, data[0,:] - data_star[0,:] - data_t[0,:])

``````
``````

In [ ]:

plt.hist(data[0,:] - data_star[0,:] - data_t[0,:])

``````
``````

In [ ]:

``````