288 lines
8.1 KiB
Python
288 lines
8.1 KiB
Python
#!/usr/bin/env python3
|
|
# -*- coding: utf-8 -*-
|
|
"""
|
|
abgewandelt von:
|
|
Created on Fri Dec 3 08:49:07 2021
|
|
|
|
@author: dreizler
|
|
"""
|
|
import pdb
|
|
|
|
from scipy.stats import gamma,norm,beta,truncnorm
|
|
import numpy as np
|
|
|
|
from types import SimpleNamespace
|
|
|
|
import dynesty
|
|
from dynesty.utils import resample_equal
|
|
|
|
# system functions that are always useful to have
|
|
import time, sys, os
|
|
|
|
# basic numeric setup
|
|
import numpy as np
|
|
|
|
from dynesty import plotting as dyplot
|
|
|
|
from oct2py import Oct2Py
|
|
|
|
from multiprocessing import Pool
|
|
|
|
# plotting
|
|
import matplotlib
|
|
from matplotlib import pyplot as plt
|
|
# inline plotting
|
|
#%matplotlib inline
|
|
|
|
import csv
|
|
import json
|
|
import pickle
|
|
|
|
def parse_parameters(settings):
|
|
p = []
|
|
fix = {}
|
|
for planet in settings["planets"]:
|
|
for key, value in planet.items():
|
|
fix[key] = [(None, None)] * len(settings["planets"])
|
|
|
|
i = 0
|
|
for planet in settings["planets"]:
|
|
for key, value in planet.items():
|
|
if value["type"] != "fixed":
|
|
p += [(i, key, value["type"], value["values"])]
|
|
else:
|
|
fix[key][i] = value["values"]
|
|
i += 1
|
|
|
|
for key, value in settings["parameters"].items():
|
|
if value["type"] != "fixed":
|
|
p += [(-1, key, value["type"], value["values"])]
|
|
else:
|
|
fix[key] = value["values"]
|
|
|
|
return p, fix
|
|
|
|
def parse_settings():
|
|
settings = json.load(open("parameter.json"))
|
|
|
|
p, fix = parse_parameters(settings)
|
|
|
|
return p, fix, settings["nkern"], settings["nlivepoints"]
|
|
|
|
|
|
def generate_mockdata():
|
|
np.random.seed(56101)
|
|
|
|
P_true=[12.3456, 29.3]
|
|
Tmid0_true=[9.8765432, 15.4029]
|
|
ror_true=[0.021, 0.027]
|
|
aor_true=30
|
|
mu_true=[1.24e-5, 3e-5]
|
|
ex0_true=[0.01234, 0.00998]
|
|
ey0_true=[0.02345, 0.0189]
|
|
I0_true=[0.01745329, 0.0203]
|
|
Omega0_true=[0.4014257, 0.369]
|
|
sig_true = np.log(0.000084)
|
|
|
|
#t_mock = np.linspace(0,100,500)
|
|
u1=0.5;
|
|
u2=0.3;
|
|
t=np.linspace(0,100,2) #0:0.002:100;
|
|
tRV=np.linspace(0,100,1000)
|
|
|
|
init_octpy()
|
|
LC,RV_o_r,Y_o_r,Z_o_r,O = octave.feval("AnalyticLC",P_true, Tmid0_true, ror_true, aor_true, mu_true, ex0_true, ey0_true, I0_true, Omega0_true, t, u1, u2,'tRV',tRV, nout=5)
|
|
|
|
#y_old = octave.AnalyticLC(P_true, Tmid0_true, ror_true, aor_true, mu_true, ex0_true, ey0_true, I0_true, Omega0_true, t, u1, u2)
|
|
y_true2 = []
|
|
for i in RV_o_r:
|
|
y_true2.append(i[0])
|
|
y_true2 = np.array(y_true2)
|
|
y_true = y_true2[(np.isfinite(y_true2))]
|
|
|
|
y = y_true + 0.0001*np.random.randn(len(y_true)) #0.0001
|
|
yerr = 0.00005 + np.zeros(len(y)) #0.00005
|
|
|
|
# plot results
|
|
plt.figure(figsize=(10, 5))
|
|
plt.errorbar(tRV, y, yerr=yerr, fmt='ko', ecolor='red')
|
|
plt.plot(tRV, y_true, color='blue', lw=3)
|
|
plt.xlabel(r'$t$')
|
|
plt.ylabel(r'$y_true$')
|
|
plt.tight_layout()
|
|
plt.show()
|
|
plt.close()
|
|
|
|
return y, yerr, y_true
|
|
|
|
def calculate(nkern, nlivepoints, parameters, fixed_parameters):
|
|
ndim = len(parameters)
|
|
p = Pool(nkern, init, [parameters, fixed_parameters])
|
|
|
|
#P, Tmid0, ror, aor, mu, ex0, ey0, I0, Omega0, sig
|
|
dsampler = dynesty.NestedSampler(loglike, prior_transform, ndim=ndim, periodic=None,
|
|
bound='multi', sample='auto',nlive=nlivepoints,pool=p, queue_size=nkern) #n*n/2 live points mind.
|
|
dsampler.run_nested(dlogz=0.001)
|
|
dres = dsampler.results
|
|
with open('data_4dim.pickle', 'wb') as f:
|
|
pickle.dump(dres, f, pickle.HIGHEST_PROTOCOL)
|
|
|
|
return dres
|
|
|
|
def create_plot(dres):
|
|
P_true=[12.3456, 29.3]
|
|
Tmid0_true=[9.8765432, 15.4029]
|
|
ror_true=[0.021, 0.027]
|
|
aor_true=30
|
|
mu_true=[1.24e-5, 3e-5]
|
|
ex0_true=[0.01234, 0.00998]
|
|
ey0_true=[0.02345, 0.0189]
|
|
I0_true=[0.01745329, 0.0203]
|
|
Omega0_true=[0.4014257, 0.369]
|
|
sig_true = np.log(0.000084)
|
|
|
|
truths = [P_true[0], P_true[1], sig_true]
|
|
labels = [r'P $d$',r'P $d$', r'T $d$',r'T $d$', r'ex0', r'ex0', r'ey0', r'ey0',r'$\sigma \log(m/s)$']
|
|
fig, axes = dyplot.traceplot(dres, truths=truths, truth_color='black', labels=labels, connect=True,
|
|
fig=plt.subplots(ndim, 2, constrained_layout=True))
|
|
fig.tight_layout()
|
|
plt.savefig("tmp_2p/plot_test1.png")
|
|
fig, axes = dyplot.cornerplot(dres, truths=truths, show_titles=True,
|
|
title_kwargs={'y': 1.04}, labels=labels,
|
|
fig=plt.subplots(ndim, ndim,constrained_layout=True))
|
|
|
|
fig.tight_layout()
|
|
plt.savefig("tmp_2p/plot_test2.png")
|
|
#lt.close()
|
|
|
|
def init(parameters, fixed_parameters):
|
|
global p
|
|
p = parameters
|
|
global fix
|
|
fix = fixed_parameters
|
|
init_octpy()
|
|
|
|
def init_octpy():
|
|
global octave
|
|
octave = Oct2Py()
|
|
BaseDir='/home/end/anal/endsversion/'
|
|
octave.addpath(BaseDir)
|
|
octave.addpath(BaseDir + 'AnalyticLC/')
|
|
octave.addpath(BaseDir + 'FittingAndMath/')
|
|
octave.addpath(BaseDir + 'FluxCalculation/')
|
|
octave.LaplaceCoeffs('load',BaseDir + 'LaplaceCoeffs.mat')
|
|
|
|
# Prior transforms for nested samplers:
|
|
def transform_uniform(x, hyperparameters):
|
|
a, b = hyperparameters
|
|
return a + (b-a)*x
|
|
|
|
def transform_loguniform(x, hyperparameters):
|
|
a, b = hyperparameters
|
|
la = np.log(a)
|
|
lb = np.log(b)
|
|
return np.exp(la + x * (lb - la))
|
|
|
|
def transform_normal(x, hyperparameters):
|
|
mu, sigma = hyperparameters
|
|
return norm.ppf(x, loc=mu, scale=sigma)
|
|
|
|
def transform_beta(x, hyperparameters):
|
|
a, b = hyperparameters
|
|
return beta.ppf(x, a, b)
|
|
|
|
def transform_exponential(x, hyperparameters):
|
|
a = hyperparameters
|
|
return gamma.ppf(x, a)
|
|
|
|
def transform_truncated_normal(x, hyperparameters):
|
|
mu, sigma, a, b = hyperparameters
|
|
ar, br = (a - mu) / sigma, (b - mu) / sigma
|
|
return truncnorm.ppf(x, ar, br, loc=mu, scale=sigma)
|
|
|
|
def transform_modifiedjeffreys(x, hyperparameters):
|
|
turn, hi = hyperparameters
|
|
return turn * (np.exp( (x + 1e-10) * np.log(hi/turn + 1)) - 1)
|
|
|
|
# prior transform
|
|
def prior_transform(utheta):
|
|
global p;
|
|
global fix;
|
|
|
|
transforms = [x[2] for x in p]
|
|
hyperparams = [x[3] for x in p]
|
|
transformed_priors = np.zeros(len(utheta))
|
|
for k,param in enumerate(utheta):
|
|
if transforms[k] == 'uniform':
|
|
transformed_priors[k] = transform_uniform(param,hyperparams[k])
|
|
elif transforms[k] == 'loguniform':
|
|
transformed_priors[k] = transform_loguniform(param,hyperparams[k])
|
|
elif transforms[k] == 'normal':
|
|
transformed_priors[k] = transform_normal(param,hyperparams[k])
|
|
else:
|
|
raise ValueError("Transform not known")
|
|
|
|
return transformed_priors
|
|
|
|
# log-likelihood
|
|
def loglike(theta):
|
|
global octave;
|
|
global p;
|
|
global fix;
|
|
parameters = fix.copy()
|
|
|
|
for i in range(len(theta)):
|
|
planetid = p[i][0]
|
|
keyname = p[i][1]
|
|
if planetid == -1:
|
|
parameters[keyname] = theta[i]
|
|
else:
|
|
parameters[keyname][planetid] = theta[i]
|
|
|
|
n = SimpleNamespace(**parameters)
|
|
|
|
u1=0.5;
|
|
u2=0.3;
|
|
t=np.linspace(0,100,2) #0:0.002:100;
|
|
tRV=np.linspace(0,100,1000)
|
|
|
|
try:
|
|
LC,RV_o_r,Y_o_r,Z_o_r,O = octave.feval("AnalyticLC",n.P, n.Tmid0, n.ror, n.aor, n.mu, n.ex0, n.ey0, n.I0, n.Omega0, t, u1, u2,'tRV',tRV, nout=5)
|
|
|
|
model =[]
|
|
for i in RV_o_r:
|
|
model.append(i[0])
|
|
model = np.array(model)
|
|
model23 = model[~(np.isnan(model))]
|
|
model2 = model23[~(np.iscomplex(model23))]
|
|
#print("Model:")
|
|
#print(model2)
|
|
#print("y:")
|
|
#print(y)
|
|
if (model2.size == 0):
|
|
#print("error")
|
|
return -np.inf
|
|
inv_sigma2 = 1.0 / (n.yerr**2 + np.exp(2 * n.sig))
|
|
ret = -0.5 * (np.sum((n.y-model2)**2 * inv_sigma2 - np.log(inv_sigma2)))
|
|
#if np.imag(ret)!=0:
|
|
# return -np.inf
|
|
#print(ret)
|
|
return ret
|
|
except Exception as e:
|
|
print(e)
|
|
return -np.inf
|
|
|
|
def main():
|
|
parameters, fixed, nkern, nlivepoints = parse_settings()
|
|
|
|
y, yerr, y_true = generate_mockdata()
|
|
|
|
fixed["y"] = y
|
|
fixed["yerr"] = yerr
|
|
fixed["y_true"] = y_true
|
|
|
|
results = calculate(nkern, nlivepoints, parameters, fixed)
|
|
create_plot(results)
|
|
|
|
main()
|