Commit 5f9441cc authored by Sweta Shah's avatar Sweta Shah
Browse files

fast fisher matrix scripts

parent 0ef23ec8
import FastBinary
from synthlisa import *
import numpy as np
import matplotlib.pyplot as plt
import sys
sys.path.append("..")
import read_catalogue
from compute_foreground_signal import runFastBinarycode
from compute_SNR_catalogue import calc_SNR_6links_AE
from foreground_utils import *
########## THE configuration parameters for the simulation
OneYear = 2**25
NYear = 3
stime = 15 ; patches = 256
samples = int(NYear*OneYear/stime)
## Lisasolve configurations
armL = 2.5e9
## Compute the Instrumental noise
c = 2.99792458e8
Spm = 3.0e-15**2/(4.0*math.pi*math.pi*c**2) ## 2.5e-48
Sop = ((10.0e-12)**2)*(4.0*math.pi*math.pi/c**2) ## 4.4e-38
print 'Spm, Sop', Spm, Sop
Sl = 0
lisa = EccentricInclined(2.5e9/2.99792458e8,0.0,0.,-1,0)
noiseTDI = TDInoise(lisa, stime,Spm, stime,Sop, stime,Sl)
[tn, Xmn, Ymn, Zmn] = np.transpose(synthlisa.getobs(samples,stime,[noiseTDI.t, noiseTDI.Xm, noiseTDI.Ym, noiseTDI.Zm]))
## Converting AlphaBetaGamma->AET for the instrumental noise as SyntheticLISA does have option for outputs alpha, beta, gamma
Amn, Emn, Tmn = xyz_to_AET(Xmn, Ymn, Zmn)
## The PSD of the instrumental noise
specAn = spect(Amn,stime,patches)
specEn = spect(Emn,stime,patches)
Inst_noise = [specAn[1:,0], specAn[1:,1], specEn[1:,1]]
SNR_threshold = 7.0
criteria_stop_loop = 100
######## THE reduced data file which is *not* a table ########
[freq, fdot, elat, elon, amp, inc, psi, phi0] = np.load("/Users/swshah/Desktop/simulations_lisasolve/my_scripts/binary_catalogue/Toonen_Aug2017/ag/ST_ag_catalogue_Full_GW_table.npy")
## Set the path for saving all the files here
savepath = "/Users/swshah/Desktop/simulations_lisasolve/my_scripts/catalogue_simulation_results/Toonen_Aug2017/ag_6links/3yr_15s/"
iter_num = 0
num_sources_SNR_gt_thres = len(amp)
while (num_sources_SNR_gt_thres > criteria_stop_loop):
## Rewrite the old values in mytable to that of the new ones
GW_binary_parameters = [freq, fdot, elat, elon, amp, inc, psi, phi0]
## Rewrite the old values in GW_binary_parameters to that of the new ones
mytable = np.transpose(GW_binary_parameters)
## Compute foreground signal
print "Compute the" + str(iter_num) + " foreground signal"
Xc,Yc,Zc, timet = runFastBinarycode(mytable, iter_num, OneYear, NYear, stime, armL, savepath)
## Converting XYZ->AET for the confusion noise because FastBinary does not output alpha, beta gamma, rather XYZ
Ac, Ec, Tc = xyz_to_AET(Xc,Yc,Zc)
## Below we combine foreground signal AET (using conversion from XYZ) and instrumental noise AET (using conversion from AlphaBetaGamma)
if (len(Amn) == len(Ac)):
specAnc = spect(Ac+Amn, stime, patches)
specEnc = spect(Ec+Emn, stime, patches)
else:
specAnc = spect(Ac+Amn[0:len(Ac)], stime, patches)
specEnc = spect(Ec+Emn[0:len(Ac)], stime, patches)
InstConf_noise = [specAnc[1:,0], specAnc[1:,1], specEnc[1:,1]]
#plt.figure(1)
#plt.subplot(221)
#plt.loglog(specXnc[1:,0], specXnc[1:,1])
#plt.show(block=False)
## Compute their SNR
SNR_catalogue_AE = calc_SNR_6links_AE(GW_binary_parameters, iter_num, armL, NYear, OneYear, stime, InstConf_noise, savepath)
print "Done computing the SNR of "+ str(len(SNR_catalogue_AE)) + " sources"
print "max SNR", round(max(SNR_catalogue_AE))
#plt.subplot(222)
#plt.hist(np.log10(SNR_catalogue), normed=True); plt.axvline(math.log10(20.0), color='black'); plt.title('Log SNR'); plt.show(block=False)
## Filter the sources with SNR < threshold
idx_lt_20 = np.where(SNR_catalogue_AE < SNR_threshold)[0]
idx_gt_20 = np.where(SNR_catalogue_AE > SNR_threshold)[0]
## Update the number of sources with SNR that are less than the threshold
num_sources_SNR_lt_thres = len(idx_lt_20)
num_sources_SNR_gt_thres = len(idx_gt_20)
#print "--- num of sources with SNR < threhold ", num_sources_SNR_lt_thres
print "--- num of sources with SNR > threhold ", num_sources_SNR_gt_thres
## Update mytable with U N R E S O L V A B L E sources
amp_lt_20 = amp[idx_lt_20]
freq_lt_20 = freq[idx_lt_20]
fdot_lt_20 = fdot[idx_lt_20]
elat_lt_20 = elat[idx_lt_20]
elon_lt_20 = elon[idx_lt_20]
inc_lt_20 = inc[idx_lt_20]
psi_lt_20 = psi[idx_lt_20]
phi0_lt_20 = phi0[idx_lt_20]
#plt.subplot(223)
#plt.hist(np.log10(amp_lt_20), normed=True); plt.title('Log Amp'); plt.show(block=False)
##Update the parameters with the new number of reduced values
GW_binary_parameters_lt_20 = [freq_lt_20, fdot_lt_20, elat_lt_20, elon_lt_20, amp_lt_20, inc_lt_20, psi_lt_20, phi0_lt_20]
mytable_lt_20 = np.transpose(GW_binary_parameters_lt_20)
#counter = counter +1
#print "A F T E R W A R D S Shape of updated table ==>", np.shape(mytable_lt_20)
#print " Finished after loop number ", iter_num
#filename_unresolved_sources = str(int(NYear))+"yr_"+str(stime)+"s_L_"+str(int(armL))[0:2]+"unresolved_sources_after_round"+str(iter_num)+".npy"
#np.save(savepath+"unresolved_binaries/"+filename_unresolved_sources, GW_binary_parameters_lt_20)
############### R E S O L V E D ###############
## get those sources that are Resolvable
amp_gt_20 = amp[idx_gt_20]
freq_gt_20 = freq[idx_gt_20]
fdot_gt_20 = fdot[idx_gt_20]
elat_gt_20 = elat[idx_gt_20]
elon_gt_20 = elon[idx_gt_20]
inc_gt_20 = inc[idx_gt_20]
psi_gt_20 = psi[idx_gt_20]
phi0_gt_20 = phi0[idx_gt_20]
#plt.subplot(224)
#plt.hist(np.log10(amp_gt_20), normed=True, histtype='step'); plt.title('Log Amp'); plt.show(block=False)
GW_binary_parameters_gt_20 = [freq_gt_20, fdot_gt_20, elat_gt_20, elon_gt_20, amp_gt_20, inc_gt_20, psi_gt_20, phi0_gt_20]
filename_resolved_sources = str(int(NYear))+"yr_"+str(stime)+"s_L_"+str(int(armL))[0:2]+"resolved_sources_after_round"+str(iter_num)+".npy"
np.save(savepath+"resolved_binaries/"+filename_resolved_sources, GW_binary_parameters_gt_20)
## Rewrite the old values in individual parameters to that of the new ones
amp = amp_lt_20
freq = freq_lt_20
fdot = fdot_lt_20
elat = elat_lt_20
elon = elon_lt_20
inc = inc_lt_20
psi = psi_lt_20
phi0 = phi0_lt_20
iter_num = iter_num +1
################# S T A T I S T I C S #################
#--- num of sources with SNR > threhold
# 3yr
#Compute the0 foreground signal
#Done computing the SNR of 22159689 sources max SNR 688.0
#--- num of sources with SNR > threhold 7235
#Compute the1 foreground signal
#Done computing the SNR of 22152454 sources max SNR 43.0
#--- num of sources with SNR > threhold 6800
#Compute the2 foreground signal
#Done computing the SNR of 22145654 sources. max SNR 28.0
#--- num of sources with SNR > threhold 3870
#Compute the3 foreground signal
#Done computing the SNR of 22141784 sources. max SNR 12.0
#--- num of sources with SNR > threhold 2052
#Compute the4 foreground signal
#Done computing the SNR of 22139732 sources. max SNR 10.0
#--- num of sources with SNR > threhold 1118
#Compute the5 foreground signal
#Done computing the SNR of 22138614 sources. max SNR 9.0
#--- num of sources with SNR > threhold 661
#Compute the6 foreground signal
#Done computing the SNR of 22137953 sources max SNR 9.0
#--- num of sources with SNR > threhold 417
#Compute the7 foreground signal
#Done computing the SNR of 22137536 sources. max SNR 8.0
#--- num of sources with SNR > threhold 258
#Compute the8 foreground signal
#Done computing the SNR of 22137278 sources. max SNR 8.0
#--- num of sources with SNR > threhold 178
#Compute the9 foreground signal
#Done computing the SNR of 22137100 sources. max SNR 8.0
#--- num of sources with SNR > threhold 139
#Compute the10 foreground signal
#Done computing the SNR of 22136961 sources. max SNR 8.0
#--- num of sources with SNR > threhold 97
sys.exit()
#--- num of sources with SNR > threhold 1808
# 7235
# 6800
# 3870
# 2052
# 1118
# 661
# 417
# 258
# 178
# 139
# 97
#[7235, 6800, 3870, 2052, 1118, 661, 417, 258, 178, 139, 97]
## 22825
import FastBinary
import matplotlib.pyplot as plt
plt.rcParams['agg.path.chunksize'] = 10000
import numpy as np
from math import *
import sys
sys.path.append("..")
from formula_utils import *
from fast_fisher_utils import *
import time
start_time = time.time()
## Detector configuration
armL = 2.5e9
FastBinary.setL(armL)
fastbin = FastBinary.FastGalacticBinary()
OneYear = 2**25
NYear = 4
stime = 15 ; patches = 256
headpath = "/Users/swshah/Desktop/simulations_lisasolve/my_scripts/catalogue_simulation_results/Toonen_Aug2017/ag_6links/4yr_15s/"
f_bins, fdots, elats, elons, ampls, incs, psis, phi0s = [], [], [], [], [], [], [], []
## Load the parameters of all the resolved binaries ##
iter_nums = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
num_res_bins = []
for i in range(len(iter_nums)):
iter_num = iter_nums[0]
filename_resolved_sources = str(int(NYear))+"yr_"+str(stime)+"s_L_"+str(int(armL))[0:2]+"resolved_sources_after_round"+str(iter_nums[i])+".npy"
GW_par_res = np.load(headpath+"resolved_binaries/"+filename_resolved_sources)
print iter_nums[i], len(GW_par_res[0])
num_res_bins.append(len(GW_par_res[0]))
f_bins.extend(GW_par_res[0]); fdots.extend(GW_par_res[1])
elats.extend(GW_par_res[2]); elons.extend(GW_par_res[3])
ampls.extend(GW_par_res[4]); incs.extend(GW_par_res[5])
psis.extend(GW_par_res[6]); phi0s.extend(GW_par_res[7])
tot_res_bins = sum(num_res_bins)
print "total number of resolved sources:", tot_res_bins ## 27547
#0 9029
#1 8222
#2 4721
#3 2504
#4 1353
#5 785
#6 429
#7 264
#8 159
#9 81
import synthlisa as synthlisa
## Noise
c = 2.99792458e8
Spm = 3.0e-15**2/(4.0*math.pi*math.pi*c**2) ## 2.5e-48
Sop = ((10.0e-12)**2)*(4.0*math.pi*math.pi/c**2) ## 4.4e-38
Sl = 0 # no laser noise since it would be canceled with first-generation TDI
samples = int(NYear*OneYear/stime)
lisa = synthlisa.EccentricInclined(2.5e9/2.99792458e8,0.0,0.,-1,0)
noiseTDI = synthlisa.TDInoise(lisa, stime,Spm, stime,Sop, stime,Sl)
[tn, Xmn, Ymn, Zmn] = np.transpose(synthlisa.getobs(samples,stime,[noiseTDI.t, noiseTDI.Xm, noiseTDI.Ym, noiseTDI.Zm]))
##### Confusion Noise - use the last instance (Reduced) #####
conf = np.load(headpath+"background_signal/"+"4yr_15s_time_L_25_Xm_Xc_after_round9.npy")
Xc, Yc, Zc, tt = conf[0], conf[1], conf[2], conf[3]
if len(Xmn) > len(Xc):
Xmn = Xmn[0:len(Xc)]
Ymn = Ymn[0:len(Yc)]
Zmn = Zmn[0:len(Zc)]
An, En, Tn = AET_xyz(Xmn,Ymn,Zmn)
specAn = spect(An,stime,patches)
specEn = spect(En,stime,patches)
plt.loglog(specAn[1:,0], specAn[1:,1], color='b', ls='--')
plt.loglog(specAn[1:,0], specEn[1:,1], color='g', ls='--')
Ac, Ec, Tc = AET_xyz(Xc,Yc,Zc)
nspec = [synthlisa.spect(An+Ac,stime,patches), synthlisa.spect(En+Ec,stime,patches)]
ep = sqrt(2.220446049250313e-16)
pre_f = 1e0*ep
delat = 1e-3
delon = 1e-3
df, de, dn = [], [], []
SNRs = []
for id_binary in range(tot_res_bins):
## Load each binary parameter
f_bin, fdot, ampl = f_bins[id_binary], fdots[id_binary], ampls[id_binary]
df_bin = pre_f*f_bin
inc, elat, elon = incs[id_binary], elats[id_binary], elons[id_binary]
phi0, psi = phi0s[id_binary], psis[id_binary]
##
fastbin.Frequency = f_bin
fastbin.FrequencyDerivative = fdot
fastbin.EclipticLatitude = elat
fastbin.EclipticLongitude = elon
fastbin.Amplitude = ampl
fastbin.Inclination = inc
fastbin.Polarization = psi
fastbin.InitialPhase = phi0
(hf,Yf,Zf) = fastbin.onefourier(T=NYear*OneYear, dt=stime)
(ht,Yt,Zt) = fastbin.TDI(T=NYear*OneYear,dt=stime)
## Convert to the optimal signals
At, Et, Tt = AET_xyz(ht,Yt,Zt)
##snrs
SNR_At = synthlisa.sn(At, An+Ac, stime, 256)
SNR_Et = synthlisa.sn(Et, En+Ec, stime, 256)
SNR_tot_temp = sqrt(SNR_At**2 + SNR_Et**2)
SNRs.append(SNR_tot_temp)
## Do the fisher analysis
temp = f_bin+df_bin; df_bin = temp-f_bin
temp5 = elat+delat; delat = temp5-elat
temp6 = elon+delon; delon = temp6-elon
##
params = [f_bin, fdot, phi0, inc, ampl, elat, elon, psi]
dxs = [df_bin, delat, delon]
fisher = my_calc_fisher_3X3_symm(params, dxs, nspec, NYear, OneYear, stime, 'time', 'A','E')
cov = cov_3X3(fisher)
df.append(cov[0]), de.append(cov[3]), dn.append(cov[5])
#print id_binary, SNR_tot_temp, cov[0]/f_bin,3, cov[3], cov[5]
print id_binary, round(SNR_tot_temp,2), round(cov[0]/f_bin,14), round(cov[3],5), round(cov[5],5)
print "Computed in ==> ", (time.time() - start_time)
# Computed in ==>
sys.exit()
savepath = "/Users/swshah/LISAGalacticBinaryDataAnalysis/fast_fisher_analysis/GBs_2pt5/results/"
np.save(savepath+"GB_fisher_analysis_freq_sky_4yr_6link", [SNRs, df, de, dn
])
sys.exit()
#!/usr/bin/env python
import FastBinary
from fast_signals_utils import *
from math import sqrt
def my_calc_fisher_2X2_symm(params, dxs, nspec, NYear, OneYear, stime, dom, *obs):
freq, fdot, phi0, inc, ampl = params[0], params[1], params[2], params[3], params[4]
elat, elon, psi = params[5], params[6], params[7]
##
dfreq, dphi, dinc, dampl = dxs[0], dxs[1], dxs[2], dxs[3]
delat, delon, dpsi = dxs[4], dxs[5], dxs[6]
## Signals
s_nA_p = sig_nA(params, ampl+dampl, NYear, OneYear, stime, dom, *obs)
s_ninc_p = sig_ninc(params, inc+dinc, NYear, OneYear, stime, dom, *obs)
#print "s_nA_p, s_ninc_p", s_nA_p.shape, s_ninc_p.shape
s_nA_m = sig_nA(params, ampl-dampl, NYear, OneYear, stime, dom, *obs)
s_ninc_m = sig_ninc(params, inc-dinc, NYear, OneYear, stime, dom, *obs)
#print "s_nA_m, s_ninc_m", s_nA_m.shape, s_ninc_m.shape
## Partial derivatives with respect to corresponding parameters
print "fast_fisher_utils, len(s_nA_m), len(s_ninc_m)", len(s_nA_m), len(s_ninc_m)
deriv_A = dh_dx_symm(s_nA_p, s_nA_m, dampl, *obs)
deriv_cos_inc = (-1./(np.sin(inc)))*np.array((dh_dx_symm(s_ninc_p, s_ninc_m, dinc, *obs)))
print "len(deriv_A), len(deriv_cos_inc)", len(deriv_A), len(deriv_cos_inc)
## fisher matrix elements
print "prep_innerprod"
gaa = prep_innerprod(deriv_A, deriv_A, nspec, stime, *obs)
gii = prep_innerprod(deriv_cos_inc, deriv_cos_inc, nspec, stime, *obs)
gai = prep_innerprod(deriv_A, deriv_cos_inc, nspec, stime, *obs)
#print gaa, gai, gii
retval = [gaa, gai, gii]
return np.array(retval)
def my_calc_fisher_2X2_symm_freq(params, dxs, nspec, NYear, OneYear, stime, df, *obs):
freq, fdot, phi0, inc, ampl = params[0], params[1], params[2], params[3], params[4]
elat, elon, psi = params[5], params[6], params[7]
##
dfreq, dphi, dinc, dampl = dxs[0], dxs[1], dxs[2], dxs[3]
delat, delon, dpsi = dxs[4], dxs[5], dxs[6]
## Signals
s_nA_p = sigF_nA(params, ampl+dampl, NYear, OneYear, stime, *obs)
s_ninc_p = sigF_ninc(params, inc+dinc, NYear, OneYear, stime, *obs)
#print "s_nA_p, s_ninc_p", s_nA_p.shape, s_ninc_p.shape
s_nA_m = sigF_nA(params, ampl-dampl, NYear, OneYear, stime, *obs)
s_ninc_m = sigF_ninc(params, inc-dinc, NYear, OneYear, stime, *obs)
#print "s_nA_m, s_ninc_m", s_nA_m.shape, s_ninc_m.shape
## Partial derivatives with respect to corresponding parameters
deriv_A = dh_dx_symm(s_nA_p, s_nA_m, dampl, *obs)
deriv_cos_inc = (-1./(np.sin(inc)))*(dh_dx_symm(s_ninc_p, s_ninc_m, dinc, *obs))
#print "deriv_A, deriv_cos_inc", deriv_A.shape, deriv_cos_inc.shape, nspec.shape
## fisher matrix elements
#print deriv_A.shape, nspec.shape
gaa = prep_innerprod_freq(deriv_A, deriv_A, nspec, stime, df, *obs)
gii = prep_innerprod_freq(deriv_cos_inc, deriv_cos_inc, nspec, stime, df, *obs)
gai = prep_innerprod_freq(deriv_A, deriv_cos_inc, nspec, stime, df, *obs)
#print gaa, gai, gii
retval = [gaa, gai, gii]
return np.array(retval)
def my_calc_fisher_7X7_symm_freq(params, dxs, nspec, NYear, OneYear, stime, dom, *obs):
freq, fdot, phi0, inc, ampl = params[0], params[1], params[2], params[3], params[4]
elat, elon, psi = params[5], params[6], params[7]
##
dfreq, dphi, dinc, dampl = dxs[0], dxs[1], dxs[2], dxs[3]
delat, delon, dpsi = dxs[4], dxs[5], dxs[6]
s_nphi_p = sig_nphi(params, phi0+dphi, NYear, OneYear, stime, dom, *obs)
s_nfreq_p = sig_nfreq(params, freq+dfreq, NYear, OneYear, stime, dom, *obs)
s_nA_p = sig_nA(params, ampl+dampl, NYear, OneYear, stime, dom, *obs)
s_ninc_p = sig_ninc(params, inc+dinc, NYear, OneYear, stime, dom, *obs)
s_npsi_p = sig_npsi(params, psi+dpsi, NYear, OneYear, stime, dom, *obs)
s_nelat_p = sig_nelat(params, elat+delat, NYear, OneYear, stime, dom, *obs)
s_nelon_p = sig_nelon(params, elon+delon, NYear, OneYear, stime, dom, *obs)
s_nphi_m = sig_nphi(params, phi0-dphi, NYear, OneYear, stime, dom, *obs)
s_nfreq_m = sig_nfreq(params, freq-dfreq, NYear, OneYear, stime, dom, *obs)
s_nA_m = sig_nA(params, ampl-dampl, NYear, OneYear, stime, dom, *obs)
s_ninc_m = sig_ninc(params, inc-dinc, NYear, OneYear, stime, dom, *obs)
s_npsi_m = sig_npsi(params, psi-dpsi, NYear, OneYear, stime, dom, *obs)
s_nelat_m = sig_nelat(params, elat-delat, NYear, OneYear, stime, dom, *obs)
s_nelon_m = sig_nelon(params, elon-delon, NYear, OneYear, stime, dom, *obs)
## Partial derivatives with respect to corresponding parameters
deriv_A = dh_dx_symm(s_nA_p, s_nA_m, dampl, *obs)
deriv_cos_inc = (-1./(np.sin(inc)))*(dh_dx_symm(s_ninc_p, s_ninc_m, dinc, *obs))
deriv_freq = dh_dx_symm(s_nfreq_p, s_nfreq_m, dfreq, *obs)
deriv_phi = dh_dx_symm(s_nphi_p, s_nphi_m, dphi, *obs)
deriv_psi = dh_dx_symm(s_npsi_p, s_npsi_m, dpsi, *obs)
deriv_sin_elat = (1./(np.cos(elat)))*(dh_dx_symm(s_nelat_p, s_nelat_m, delat, *obs))
deriv_elon = dh_dx_symm(s_nelon_p, s_nelon_m, delon, *obs)
## fisher matrix elements
gaa = prep_innerprod_freq(deriv_A, deriv_A, nspec, stime, *obs)
gii = prep_innerprod_freq(deriv_cos_inc, deriv_cos_inc, nspec, stime, *obs)
gff = prep_innerprod_freq(deriv_freq, deriv_freq, nspec, stime, *obs)
gpp = prep_innerprod_freq(deriv_phi, deriv_phi, nspec, stime, *obs)
gss = prep_innerprod_freq(deriv_psi, deriv_psi, nspec, stime, *obs)
gee = prep_innerprod_freq(deriv_sin_elat, deriv_sin_elat, nspec, stime, *obs)
gnn = prep_innerprod_freq(deriv_elon, deriv_elon, nspec, stime, *obs)
# MM = [gaa, gap, gai, gaf, gas, gae, gan
# gpa, gpp, gpi, gpf, gps, gpe, gpn
# gia, gip, gii, gif, gis, gie, gin
# gfa, gfp, gfi, gff, gfs, gfe, gfn
# gsa, gsp, gsi, gsf, gss, gse, gsn
# gea, gep, gei, gef, ges, gee, gen
# gna, gnp, gni, gnf, gns, gne, gnn]
## OFF diagonals
gap = prep_innerprod_freq(deriv_A, deriv_phi, nspec, stime, *obs)
gai = prep_innerprod_freq(deriv_A, deriv_cos_inc, nspec, stime, *obs)
gaf = prep_innerprod_freq(deriv_A, deriv_freq, nspec, stime, *obs)
gas = prep_innerprod_freq(deriv_A, deriv_psi, nspec, stime, *obs)
gae = prep_innerprod_freq(deriv_A, deriv_sin_elat, nspec, stime, *obs)
gan = prep_innerprod_freq(deriv_A, deriv_elon, nspec, stime, *obs)
gpi = prep_innerprod_freq(deriv_phi, deriv_cos_inc, nspec, stime, *obs)
gpf = prep_innerprod_freq(deriv_phi, deriv_freq, nspec, stime, *obs)
gps = prep_innerprod_freq(deriv_phi, deriv_psi, nspec, stime, *obs)
gpe = prep_innerprod_freq(deriv_phi, deriv_sin_elat, nspec, stime, *obs)
gpn = prep_innerprod_freq(deriv_phi, deriv_elon, nspec, stime, *obs)
gif = prep_innerprod_freq(deriv_cos_inc, deriv_freq, nspec, stime, *obs)
gis = prep_innerprod_freq(deriv_cos_inc, deriv_psi, nspec, stime, *obs)
gie = prep_innerprod_freq(deriv_cos_inc, deriv_sin_elat, nspec, stime, *obs)
gin = prep_innerprod_freq(deriv_cos_inc, deriv_elon, nspec, stime, *obs)
gfs = prep_innerprod_freq(deriv_freq, deriv_psi, nspec, stime, *obs)
gfe = prep_innerprod_freq(deriv_freq, deriv_sin_elat, nspec, stime, *obs)
gfn = prep_innerprod_freq(deriv_freq, deriv_elon, nspec, stime, *obs)
gse = prep_innerprod_freq(deriv_psi, deriv_sin_elat, nspec, stime, *obs)
gsn = prep_innerprod_freq(deriv_psi, deriv_elon, nspec, stime, *obs)
gen = prep_innerprod_freq(deriv_sin_elat, deriv_elon, nspec, stime, *obs)
retval = [gaa, gap, gai, gaf, gas, gae, gan, gpp, gpi, gpf, gps, gpe, gpn, gii, gif, gis, gie, gin, gff, gfs, gfe, gfn, gss, gse, gsn, gee, gen, gnn]
return np.array(retval)
def my_calc_fisher_7X7_symm(params, dxs, nspec, NYear, OneYear, stime, dom, *obs):
freq, fdot, phi0, inc, ampl = params[0], params[1], params[2], params[3], params[4]
elat, elon, psi = params[5], params[6], params[7]
##
dfreq, dphi, dinc, dampl = dxs[0], dxs[1], dxs[2], dxs[3]
delat, delon, dpsi = dxs[4], dxs[5], dxs[6]
## Signals
s_nphi_p = sig_nphi(params, phi0+dphi, NYear, OneYear, stime, dom, *obs)
s_nfreq_p = sig_nfreq(params, freq+dfreq, NYear, OneYear, stime, dom, *obs)
s_nA_p = sig_nA(params, ampl+dampl, NYear, OneYear, stime, dom, *obs)
s_ninc_p = sig_ninc(params, inc+dinc, NYear, OneYear, stime, dom, *obs)
s_npsi_p = sig_npsi(params, psi+dpsi, NYear, OneYear, stime, dom, *obs)
s_nelat_p = sig_nelat(params, elat+delat, NYear, OneYear, stime, dom, *obs)
s_nelon_p = sig_nelon(params, elon+delon, NYear, OneYear, stime, dom, *obs)
s_nphi_m = sig_nphi(params, phi0-dphi, NYear, OneYear, stime, dom, *obs)
s_nfreq_m = sig_nfreq(params, freq-dfreq, NYear, OneYear, stime, dom, *obs)
s_nA_m = sig_nA(params, ampl-dampl, NYear, OneYear, stime, dom, *obs)
s_ninc_m = sig_ninc(params, inc-dinc, NYear, OneYear, stime, dom, *obs)
s_npsi_m = sig_npsi(params, psi-dpsi, NYear, OneYear, stime, dom, *obs)
s_nelat_m = sig_nelat(params, elat-delat, NYear, OneYear, stime, dom, *obs)
s_nelon_m = sig_nelon(params, elon-delon, NYear, OneYear, stime, dom, *obs)
## Partial derivatives with respect to corresponding parameters
deriv_A = dh_dx_symm(s_nA_p, s_nA_m, dampl, *obs)
deriv_cos_inc = (-1./(np.sin(inc)))*np.array((dh_dx_symm(s_ninc_p, s_ninc_m, dinc, *obs)))
deriv_freq = dh_dx_symm(s_nfreq_p, s_nfreq_m, dfreq, *obs)
deriv_phi = dh_dx_symm(s_nphi_p, s_nphi_m, dphi, *obs)
deriv_psi = dh_dx_symm(s_npsi_p, s_npsi_m, dpsi, *obs)
deriv_sin_elat = (1./(np.cos(elat)))*np.array((dh_dx_symm(s_nelat_p, s_nelat_m, delat, *obs)))
deriv_elon = dh_dx_symm(s_nelon_p, s_nelon_m, delon, *obs)
## fisher matrix elements
gaa = prep_innerprod(deriv_A, deriv_A, nspec, stime, *obs)
gii = prep_innerprod(deriv_cos_inc, deriv_cos_inc, nspec, stime, *obs)
gff = prep_innerprod(deriv_freq, deriv_freq, nspec, stime, *obs)
gpp = prep_innerprod(deriv_phi, deriv_phi, nspec, stime, *obs)
gss = prep_innerprod(deriv_psi, deriv_psi, nspec, stime, *obs)
gee = prep_innerprod(deriv_sin_elat, deriv_sin_elat, nspec, stime, *obs)
gnn = prep_innerprod(deriv_elon, deriv_elon, nspec, stime, *obs)
# MM = [gaa, gap, gai, gaf, gas, gae, gan
# gpa, gpp, gpi, gpf, gps, gpe, gpn
# gia, gip, gii, gif, gis, gie, gin
# gfa, gfp, gfi, gff, gfs, gfe, gfn
# gsa, gsp, gsi, gsf, gss, gse, gsn
# gea, gep, gei, gef, ges, gee, gen
# gna, gnp, gni, gnf, gns, gne, gnn]
## OFF diagonals
gap = prep_innerprod(deriv_A, deriv_phi, nspec, stime, *obs)
gai = prep_innerprod(deriv_A, deriv_cos_inc, nspec, stime, *obs)
gaf = prep_innerprod(deriv_A, deriv_freq, nspec, stime, *obs)
gas = prep_innerprod(deriv_A, deriv_psi, nspec, stime, *obs)
gae = prep_innerprod(deriv_A, deriv_sin_elat, nspec, stime, *obs)
gan = prep_innerprod(deriv_A, deriv_elon, nspec, stime, *obs)
gpi = prep_innerprod(deriv_phi, deriv_cos_inc, nspec, stime, *obs)
</