"""Proper finite-prior RV/activity posterior, independent of the browser. Each Sobol node has equal PRIOR mass. It defines a finite hyperparameter prior, not a guarantee of converged continuous integration. All conditional linear coefficient distributions and their full covariances are retained. MIT. """ import csv from pathlib import Path import numpy as np from scipy.special import ndtri, ndtr, logsumexp from scipy.stats import qmc from rv_kernel import activity_covariance, linear_gaussian_posterior HYPER_NAMES = ["evolution_days", "rotation_days", "periodic_scale", "rv_coupling_m_s", "rv_derivative_m_s_day", "calcium_amplitude_scaled", "rv_scatter_m_s", "calcium_scatter_scaled", "second_period_days"] COEFFICIENT_NAMES = ["rv_offset_m_s", "calcium_offset_scaled", "b_sine_m_s", "b_cosine_m_s", "c_sine_m_s", "c_cosine_m_s"] def read_observations(folder): with (Path(folder) / "harps2012-corot7.csv").open() as source: rows = list(csv.DictReader(source)) arrays = {name: np.array([float(row[name]) for row in rows]) for name in rows[0]} time_origin = float(arrays["bjd_utc_day"][0]) rv_origin = 31180.0 # fixed declared reference, not a fit-dependent prior center calcium_origin = -4.72 # fixed declared reference in dimensionless log index times = arrays["bjd_utc_day"] - time_origin values = np.r_[arrays["rv_km_s"] * 1000 - rv_origin, (arrays["log_rhk"] - calcium_origin) / .01] errors = np.r_[arrays["rv_error_km_s"] * 1000, arrays["log_rhk_error"] / .01] return rows, times, values, errors, dict(time_origin_bjd_utc=time_origin, rv_origin_m_s=rv_origin, calcium_origin_log_rhk=calcium_origin, calcium_normalization_dex=.01) def prior_nodes(power, seed=1801, rotation_prior="broad", amplitude_scale=10.): uniform = qmc.Sobol(9, scramble=True, seed=seed).random_base2(power) bounded_normal = ndtri(np.clip(uniform,1e-12,1-1e-12)) log_uniform = lambda column, low, high: low * (high / low) ** uniform[:,column] rotation = 18 + 17 * uniform[:,1] if rotation_prior == "photometry": rotation = 23.81 + .03 * bounded_normal[:,1] elif rotation_prior == "wide": rotation = 3 + 57 * uniform[:,1] nodes = np.column_stack([ log_uniform(0,5.,80.), rotation, log_uniform(2,.15,2.), amplitude_scale * bounded_normal[:,3], 5 * amplitude_scale * bounded_normal[:,4], log_uniform(5,.1,3.), log_uniform(6,.3,8.), log_uniform(7,.05,2.), 3.3 + .8 * uniform[:,8]]) return nodes def design_matrix(times, second_period, planets): count = len(times) design = np.zeros((2*count,6 if planets else 2)) design[:count,0] = 1 design[count:,1] = 1 if planets: for offset, period in [(2,.85359165),(4,second_period)]: design[:count,offset] = np.sin(2*np.pi*times/period) design[:count,offset+1] = np.cos(2*np.pi*times/period) return design def posterior_at_nodes(times, values, errors, nodes, model, planet_scale=10., log_prior_minus_proposal=None): count = len(times) planets = model != "activity" dimensions = 6 if planets else 2 means = np.empty((len(nodes),dimensions)); covariances = np.empty((len(nodes),dimensions,dimensions)) log_densities = np.empty(len(nodes)); failures = [] scales = [100.,10.] + ([planet_scale]*4 if planets else []) for index, node in enumerate(nodes): parameters = dict(zip(HYPER_NAMES,node)) noise = np.diag(errors ** 2 + np.r_[np.full(count,node[6]**2),np.full(count,node[7]**2)]) if model != "kepler": noise += activity_covariance(times,times,parameters) design = design_matrix(times,node[8],planets) try: density, mean, covariance = linear_gaussian_posterior(values,noise,design,scales) if not np.all(np.isfinite(covariance)) or np.min(np.linalg.eigvalsh(covariance)) < -1e-7: raise ValueError("invalid coefficient covariance") except (ValueError,np.linalg.LinAlgError) as error: failures.append(dict(node=index,reason=str(error))) log_densities[index]=-np.inf; means[index]=0; covariances[index]=0 continue log_densities[index]=density; means[index]=mean; covariances[index]=covariance weighted_densities=log_densities if log_prior_minus_proposal is None else log_densities+log_prior_minus_proposal log_normalization=logsumexp(weighted_densities) weights=np.exp(weighted_densities-log_normalization) mean=np.einsum("n,nj->j",weights,means) centered=means-mean covariance=np.einsum("n,nij->ij",weights,covariances) + np.einsum("n,ni,nj->ij",weights,centered,centered) hyper_mean=np.einsum("n,nj->j",weights,nodes) complete_mean=np.r_[hyper_mean,mean] combined=np.column_stack([nodes,means]); centered_combined=combined-complete_mean complete_covariance=np.einsum("n,ni,nj->ij",weights,centered_combined,centered_combined) complete_covariance[9:,9:]+=np.einsum("n,nij->ij",weights,covariances) return dict(model=model,weights=weights,log_densities=log_densities, coefficient_means=means,coefficient_covariances=covariances, coefficient_mean=mean,coefficient_covariance=covariance, joint_mean=complete_mean,joint_covariance=complete_covariance, finite_prior_log_marginal=float(log_normalization-np.log(len(nodes))), effective_nodes=float(1/np.sum(weights**2)),maximum_weight=float(weights.max()),failures=failures) def conditional_draws(posterior,nodes,seed=1802,count=8192): random=np.random.default_rng(seed) chosen=random.choice(len(nodes),size=count,p=posterior["weights"]) coefficients=[] for index in chosen: covariance=posterior["coefficient_covariances"][index] eigenvalues,eigenvectors=np.linalg.eigh(covariance) if eigenvalues.min() < -1e-7: raise ValueError("posterior draw covariance is not positive") coefficients.append(posterior["coefficient_means"][index] + np.einsum("ij,j->i",eigenvectors,np.sqrt(np.maximum(0,eigenvalues))*random.normal(size=len(eigenvalues)))) return nodes[chosen],np.array(coefficients) def summarize(posterior,nodes): hyper,coefficients=conditional_draws(posterior,nodes) summary=dict(model=posterior["model"],effective_nodes=posterior["effective_nodes"], maximum_node_weight=posterior["maximum_weight"],finite_prior_log_marginal=posterior["finite_prior_log_marginal"], coefficient_names=COEFFICIENT_NAMES[:coefficients.shape[1]], coefficient_mean=posterior["coefficient_mean"].tolist(), coefficient_covariance=posterior["coefficient_covariance"].tolist(), joint_parameter_names=HYPER_NAMES+COEFFICIENT_NAMES[:coefficients.shape[1]], joint_parameter_mean=posterior["joint_mean"].tolist(),joint_parameter_covariance=posterior["joint_covariance"].tolist(), failures=posterior["failures"]) if coefficients.shape[1]==6: amplitude_b=np.hypot(coefficients[:,2],coefficients[:,3]); amplitude_c=np.hypot(coefficients[:,4],coefficients[:,5]) summary["amplitude_b_m_s_quantiles_16_50_84"]=np.quantile(amplitude_b,[.16,.5,.84]).tolist() summary["amplitude_c_m_s_quantiles_16_50_84"]=np.quantile(amplitude_c,[.16,.5,.84]).tolist() summary["period_c_day_quantiles_16_50_84"]=np.quantile(hyper[:,8],[.16,.5,.84]).tolist() return summary def transform_prior_coordinates(coordinates,rotation_prior="broad",amplitude_scale=10.): """Independent standard-normal coordinates map to the declared proper priors.""" coordinates=np.atleast_2d(coordinates) uniform=ndtr(coordinates) logarithmic=lambda index,low,high:low*(high/low)**uniform[:,index] rotation=18+17*uniform[:,1] if rotation_prior=="photometry":rotation=23.81+.03*coordinates[:,1] elif rotation_prior=="wide":rotation=3+57*uniform[:,1] return np.column_stack([logarithmic(0,5.,80.),rotation,logarithmic(2,.15,2.), amplitude_scale*coordinates[:,3],5*amplitude_scale*coordinates[:,4], logarithmic(5,.1,3.),logarithmic(6,.3,8.),logarithmic(7,.05,2.),3.3+.8*uniform[:,8]])