<jupyter_start><jupyter_text>Here's where the TensorFlow stuff starts. The "session" is the thing that will actually do that calculations.<jupyter_code>sess = tf.InteractiveSession()<jupyter_output><empty_output><jupyter_text>The way that TensorFlow works is that you define the operations on `Variable` and `placeholder` objects. Here we describe the model:<jupyter_code>T = tf.float64

# First the variables that we might want to optimize:
porb_tensor = tf.Variable(porb, dtype=T)
tp_tensor = tf.Variable(tp, dtype=T)
nu_tensor = tf.Variable(nu_arr, dtype=T)
e_param_tensor = tf.Variable(e_param, dtype=T)  # This forces the ecc to be between 0 and 1
e_tensor = 1.0 / (1.0 + tf.exp(-e_param_tensor))
varpi_tensor = tf.Variable(varpi, dtype=T)
log_sigma2_tensor = tf.Variable(0.0, dtype=T)  # Variance from observational uncertainties and model misspecification

ad_tensor = tf.Variable(a1d + np.zeros_like(nu_arr), dtype=T)

# These are some placeholders for the data:
times_tensor = tf.placeholder(T, times.shape)
dmmags_tensor = tf.placeholder(T, dmmags.shape)

# Solve Kepler's equation
mean_anom = 2.0 * np.pi * (times_tensor - tp_tensor) / porb_tensor
ecc_anom = kepler(mean_anom, e_tensor)
true_anom = 2.0 * tf.atan2(tf.sqrt(1.0+e_tensor)*tf.tan(0.5*ecc_anom), tf.sqrt(1.0-e_tensor) + tf.zeros_like(times_tensor))
# Here we define how the time delay will be calculated:
tau_tensor = -(1.0 - tf.square(e_tensor)) * tf.sin(true_anom + varpi_tensor) / (1.0 + e_tensor*tf.cos(true_anom))

# And the design matrix:
arg_tensor = 2.0 * np.pi * nu_tensor[None, :] * (times_tensor[:, None] - ad_tensor[None, :] * tau_tensor[:, None])
D_tensor = tf.concat([tf.cos(arg_tensor), tf.sin(arg_tensor)], axis=1)

# Define the linear solve for W_hat:
DTD_tensor = tf.matmul(D_tensor, D_tensor, transpose_a=True)
DTy_tensor = tf.matmul(D_tensor, dmmags_tensor[:, None], transpose_a=True)
W_hat_tensor = tf.linalg.solve(DTD_tensor, DTy_tensor)

# Finally, the model and the chi^2 objective:
model_tensor = tf.squeeze(tf.matmul(D_tensor, W_hat_tensor))
chi2_tensor = tf.reduce_sum(tf.square(dmmags_tensor - model_tensor)) * tf.exp(-log_sigma2_tensor)
chi2_tensor += len(times) * log_sigma2_tensor<jupyter_output><empty_output><jupyter_text>Here's how you could evaluate different parts of the model:<jupyter_code># We need to initialize the variables:
tf.global_variables_initializer().run()

# We'll also need to pass in the data:
data = {times_tensor: times, dmmags_tensor: dmmags}

# Let's plot the initial time delay
initial_tau = sess.run(tau_tensor, feed_dict=data)
plt.plot(times+tmid, initial_tau, ".", ms=2)
plt.ylabel(r"$\tau(t)$")
plt.xlabel("$t$");
initial_model = sess.run(model_tensor, feed_dict=data)
plt.plot(times, dmmags, ".k")
plt.plot(times, initial_model)
# plt.xlim(100, 102)
# plt.ylim(-75, 75)
plt.xlabel("t")
plt.ylabel("L(t)");<jupyter_output><empty_output><jupyter_text>Now we'll fit the parameters. We'll iterate with different subsets a few times.<jupyter_code>old_chi2 = sess.run(chi2_tensor, feed_dict=data)
for i in range(5):
    params = [log_sigma2_tensor, porb_tensor, tp_tensor]
    opt = tf.contrib.opt.ScipyOptimizerInterface(chi2_tensor, params, method="L-BFGS-B")
    opt.minimize(sess, feed_dict=data)
    
    params.append(ad_tensor)
    opt = tf.contrib.opt.ScipyOptimizerInterface(chi2_tensor, params, method="L-BFGS-B")
    opt.minimize(sess, feed_dict=data)

    params += [e_param_tensor, varpi_tensor]
    opt = tf.contrib.opt.ScipyOptimizerInterface(chi2_tensor, params, method="L-BFGS-B")
    opt.minimize(sess, feed_dict=data)
        
    new_chi2 = sess.run(chi2_tensor, feed_dict=data)
    print(old_chi2 - new_chi2)
    if np.abs(old_chi2 - new_chi2) < 1.0:
        break
    old_chi2 = new_chi2<jupyter_output>INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 38819.131340
  Number of iterations: 5
  Number of functions evaluations: 9
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 38818.046748
  Number of iterations: 9
  Number of functions evaluations: 13
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 38818.046748
  Number of iterations: 1
  Number of functions evaluations: 4
4715.05510637
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 38817.963113
  Number of iterations: 5
  Number of functions evaluations: 9
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 3[...]<jupyter_text>Let's make the updated final plots.<jupyter_code>final_tau = sess.run(tau_tensor, feed_dict=data)
plt.plot(times+tmid, initial_tau, ".", ms=2)
plt.plot(times+tmid, final_tau, ".", ms=2)
plt.ylabel(r"$\tau(t) / a$")
plt.xlabel("$t$");
models = tau_tensor[:, None] * ad_tensor[None, :]
plt.plot(times+tmid, sess.run(models, feed_dict=data), ".", ms=3)
plt.ylabel(r"$\tau(t)$")
plt.xlabel("$t$");
ivar = -np.diag(sess.run(tf.hessians(-0.5*chi2_tensor, ad_tensor), feed_dict=data)[0])
ad = sess.run(ad_tensor)
ad *= np.sign(ad[0])
sig = 1.0 / np.sqrt(ivar)
plt.errorbar(np.arange(len(ad)), ad, yerr=sig, fmt=".");
m = np.ones_like(ad, dtype=bool)
while True:
    var = 1.0 / np.sum(ivar[m])
    mu = np.sum(ivar[m] * ad[m]) * var
    var2 = np.sum(ivar[m] * (mu - ad[m])**2) * var
    m_new = np.abs(ad - mu) / np.sqrt(var) < 7.0
    if m.sum() == m_new.sum():
        m = m_new
        break
    m = m_new
ad = ad[m]
ivar = ivar[m]
m
sig = 1.0 / np.sqrt(ivar)
plt.errorbar(np.arange(len(ad)), ad, yerr=sig, fmt=".");
ad * np.sqrt(ivar)
if np.any(ad * np.sqrt(ivar) < -1.0):
    m2 = ad * np.sqrt(ivar) < -1.0
    m1 = ~m2
    ad = [
        np.sum(ivar[m1]*ad[m1]) / np.sum(ivar[m1]),
        np.sum(ivar[m2]*ad[m2]) / np.sum(ivar[m2]),
    ]
else:
    ad = [np.sum(ivar*ad) / np.sum(ivar)]
ad
if len(ad) > 1:
    inds = tf.cast(0.5 - 0.5 * (ad_tensor / tf.abs(ad_tensor)), tf.int32)
else:
    inds = tf.zeros_like(ad_tensor, dtype=tf.int32)
ad_params = tf.Variable(ad, dtype=T)
sess.run(ad_params.initializer)
ad_tensor = tf.gather(ad_params, inds)

# And the design matrix:
arg_tensor = 2.0 * np.pi * nu_tensor[None, :] * (times_tensor[:, None] - ad_tensor[None, :] * tau_tensor[:, None])
D_tensor = tf.concat([tf.cos(arg_tensor), tf.sin(arg_tensor)], axis=1)

# Define the linear solve for W_hat:
DTD_tensor = tf.matmul(D_tensor, D_tensor, transpose_a=True)
DTy_tensor = tf.matmul(D_tensor, dmmags_tensor[:, None], transpose_a=True)
W_hat_tensor = tf.linalg.solve(DTD_tensor, DTy_tensor)

# Finally, the model and the chi^2 objective:
model_tensor = tf.squeeze(tf.matmul(D_tensor, W_hat_tensor))
chi2_tensor = tf.reduce_sum(tf.square(dmmags_tensor - model_tensor)) * tf.exp(-log_sigma2_tensor)
chi2_tensor += len(times) * log_sigma2_tensor
old_chi2 = sess.run(chi2_tensor, feed_dict=data)
for i in range(5):
    params = [log_sigma2_tensor, porb_tensor, tp_tensor]
    opt = tf.contrib.opt.ScipyOptimizerInterface(chi2_tensor, params, method="L-BFGS-B")
    opt.minimize(sess, feed_dict=data)
    
    params.append(ad_params)
    opt = tf.contrib.opt.ScipyOptimizerInterface(chi2_tensor, params, method="L-BFGS-B")
    opt.minimize(sess, feed_dict=data)

    params += [e_param_tensor, varpi_tensor]
    opt = tf.contrib.opt.ScipyOptimizerInterface(chi2_tensor, params, method="L-BFGS-B")
    opt.minimize(sess, feed_dict=data)
    
    params.append(nu_tensor)
    opt = tf.contrib.opt.ScipyOptimizerInterface(chi2_tensor, params, method="L-BFGS-B")
    opt.minimize(sess, feed_dict=data)
    
    new_chi2 = sess.run(chi2_tensor, feed_dict=data)
    print(old_chi2 - new_chi2)
    if np.abs(old_chi2 - new_chi2) < 1.0:
        break
    old_chi2 = new_chi2
models = tau_tensor[:, None] * ad_tensor[None, :]
plt.plot(times+tmid, 86400.0 * sess.run(models, feed_dict=data), ".", ms=2);
plt.ylabel(r"$\tau(t)$")
plt.xlabel("$t$");
sess.run(e_tensor), e
hess_tensor = tf.hessians(-0.5*chi2_tensor, params[:-1])
hess = sess.run(hess_tensor, feed_dict=data)
1. / np.sqrt(-hess[1])
hess
np.sqrt(-np.diag(np.linalg.inv(hess[3])))
sess.run(ad_tensor)
a1d / np.sqrt(-np.diag(np.linalg.inv(hess[3])))
porb
sess.run(porb_tensor)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup