From 89b0bb282660634ecb8ce1f30e07a507d5b0a5bc Mon Sep 17 00:00:00 2001 From: Francisco Jimenez Forteza <francisco.jimenez@condor1.atlas.local> Date: Thu, 6 May 2021 16:21:09 +0000 Subject: [PATCH] rdown pe --- code_new/rdown_pe.py | 191 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 191 insertions(+) create mode 100644 code_new/rdown_pe.py diff --git a/code_new/rdown_pe.py b/code_new/rdown_pe.py new file mode 100644 index 0000000..c92e073 --- /dev/null +++ b/code_new/rdown_pe.py @@ -0,0 +1,191 @@ +# Copyright (C) 2021 Xisco Jimenez Forteza +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + +# +# ============================================================================= +# +# Preamble +# +# ============================================================================= +# +# Module to run PE on RD data +import random +from multiprocessing import Pool +import dynesty +import numpy as np +import rdown + +class Ringdown_PE: + def __init__(self,rdown_fun,data,dim,priors,errors2=1,theta=[],model='w-tau',norm_factor=0,l_int=0): + self.dim = dim + self.rdown_fun = rdown_fun + self.times = data[:,0] + self.datare = data[:,1].real + self.dataim = data[:,1].imag + self.priors = priors + self.priors_min = priors[:,0] + self.priors_max = priors[:,1] + self.prior_dim = len(priors) + self.errors2 = errors2 + self.norm_factor = norm_factor + self.model = model + self.l_int = l_int + self.theta = theta + self.dict = {'w-tau':rdown_fun.rd_model_wtau , 'w-q': rdown_fun.rd_model_wq, 'w-tau-fixed':rdown_fun.rd_model_wtau_fixed,'w-tau-fixed-m-af': rdown_fun.rd_model_wtau_m_af} + + #def log_likelihood(self,theta,sigma=1): + # """chi2 likelihood. + # """ + # modelev = dict[model](theta) + # result = -np.sum(((gwdatanew_re_tsh - modelev.real)**2+(gwdatanew_im_tsh - modelev.imag)**2)/(2*theta[-1]*error_final)) + # if np.isnan(result): + # return -np.inf + # return result + + def log_likelihood(self,theta,sigma=1): + """chi2 likelihood. + """ + modelev = self.dict[self.model](theta) + modelevre= modelev.real + modelevim= modelev.imag + + sigma2 = self.errors2 + self.l_int*(self.datare** 2+self.dataim**2) * np.exp(2 * theta[-1]) + + result = -0.5*np.sum(((self.datare - modelevre)**2+(self.dataim - modelevim)**2)/sigma2+self.l_int*(np.log(2*np.pi*sigma2)))-self.l_int*self.norm_factor + + if np.isnan(result): + return -np.inf + return result + + + def prior_transform(self,cube): + """RD uniform priors. The values for priors_min and priors_max must be given out of this function. + """ + for i in range(self.prior_dim): + cube[i] = self.priors_min[i]+ cube[i]*(self.priors_max[i]-self.priors_min[i]) + return cube + +def load_priors(model,config_parser,nmax,fitnoise=True): + # loading priors + if model == 'w-q': + tau_var_str='q' + else: + tau_var_str='tau' + + if model == 'w-tau': + w_mins=np.empty(nmax+1) + w_maxs=np.empty(nmax+1) + tau_mins=np.empty(nmax+1) + tau_maxs=np.empty(nmax+1) + a_mins=np.empty(nmax+1) + a_maxs=np.empty(nmax+1) + ph_mins=np.empty(nmax+1) + ph_maxs=np.empty(nmax+1) + + for i in range(nmax+1): + wp_min=config_parser.get('prior-w'+str(i),'w'+str(i)+'_min') + w_mins[i] = np.float(wp_min) + + wp_max=config_parser.get('prior-w'+str(i),'w'+str(i)+'_max') + w_maxs[i] = np.float(wp_max) + + taup_min=config_parser.get('prior-'+tau_var_str+str(i),tau_var_str+str(i)+'_min') + tau_mins[i] = np.float(taup_min) + + taup_max=config_parser.get('prior-'+tau_var_str+str(i),tau_var_str+str(i)+'_max') + tau_maxs[i] = np.float(taup_max) + + amp0_min=config_parser.get('prior-amp'+str(i),'amp'+str(i)+'_min') + a_mins[i] = np.float(amp0_min) + + amp1_max=config_parser.get('prior-amp'+str(i),'amp'+str(i)+'_max') + a_maxs[i] = np.float(amp1_max) + + phase_min=config_parser.get('prior-phase'+str(i),'phase'+str(i)+'_min') + ph_mins[i] = np.float(phase_min)*2*np.pi + + phase_max=config_parser.get('prior-phase'+str(i),'phase'+str(i)+'_max') + ph_maxs[i] = np.float(phase_max)*2*np.pi + + priors_min = np.concatenate((w_mins,tau_mins,a_mins,ph_mins)) + priors_max = np.concatenate((w_maxs,tau_maxs,a_maxs,ph_maxs)) + prior_dim = len(priors_min) + priors=np.column_stack((priors_min,priors_max)) + + if model == 'w-tau-fixed': + a_mins=np.empty(nmax+1) + a_maxs=np.empty(nmax+1) + ph_mins=np.empty(nmax+1) + ph_maxs=np.empty(nmax+1) + + for i in range(nmax+1): + amp0_min=config_parser.get('prior-amp'+str(i),'amp'+str(i)+'_min') + a_mins[i] = np.float(amp0_min) + + amp1_max=config_parser.get('prior-amp'+str(i),'amp'+str(i)+'_max') + a_maxs[i] = np.float(amp1_max) + + phase_min=config_parser.get('prior-phase'+str(i),'phase'+str(i)+'_min') + ph_mins[i] = np.float(phase_min)*2*np.pi + + phase_max=config_parser.get('prior-phase'+str(i),'phase'+str(i)+'_max') + ph_maxs[i] = np.float(phase_max)*2*np.pi + + + priors_min = np.concatenate((a_mins,ph_mins)) + priors_max = np.concatenate((a_maxs,ph_maxs)) + prior_dim = len(priors_min) + priors=np.column_stack((priors_min,priors_max)) + + elif model == 'w-tau-fixed-m-af': + a_mins=np.empty(nmax+1) + a_maxs=np.empty(nmax+1) + ph_mins=np.empty(nmax+1) + ph_maxs=np.empty(nmax+1) + + for i in range(nmax+1): + amp0_min=config_parser.get('prior-amp'+str(i),'amp'+str(i)+'_min') + a_mins[i] = np.float(amp0_min) + + amp1_max=config_parser.get('prior-amp'+str(i),'amp'+str(i)+'_max') + a_maxs[i] = np.float(amp1_max) + + phase_min=config_parser.get('prior-phase'+str(i),'phase'+str(i)+'_min') + ph_mins[i] = np.float(phase_min)*2*np.pi + + phase_max=config_parser.get('prior-phase'+str(i),'phase'+str(i)+'_max') + ph_maxs[i] = np.float(phase_max)*2*np.pi + + mass_min=[np.float(config_parser.get('prior-mass','mass_min'))] + mass_max=[np.float(config_parser.get('prior-mass','mass_max'))] + spin_min=[np.float(config_parser.get('prior-spin','spin_min'))] + spin_max=[np.float(config_parser.get('prior-spin','spin_max'))] + priors_min = np.concatenate((a_mins,ph_mins,mass_min,spin_min)) + priors_max = np.concatenate((a_maxs,ph_maxs,mass_max,spin_max)) + prior_dim = len(priors_min) + priors=np.column_stack((priors_min,priors_max)) + + + if fitnoise: + priors_fit_min=[np.float(config_parser.get('prior-noise','noise_min'))] + priors_fit_max=[np.float(config_parser.get('prior-noise','noise_max'))] + priors_min = np.concatenate((priors_min,priors_fit_min)) + priors_max = np.concatenate((priors_max,priors_fit_max)) + priors=np.column_stack((priors_min,priors_max)) + prior_dim = len(priors_min) + + return priors \ No newline at end of file -- GitLab