Skip to content

Instantly share code, notes, and snippets.

@evanbiederstedt
Created July 15, 2015 21:47
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save evanbiederstedt/72b0035ca4b9fce830c3 to your computer and use it in GitHub Desktop.
Save evanbiederstedt/72b0035ca4b9fce830c3 to your computer and use it in GitHub Desktop.
m^T C^-1 m values and scaling Civ
FIRST, FIND NO PEAK IN LF, NOISE PARAMETERS E-21 TO E-23
CODE:
""""
vary_x_samples125 = np.logspace(-8, -12, num=40) #num = 40
sigma125 = np.logspace(-21, -23, num=40)
Sij = vary_x_samples125[:, None, None] * norm_matrix[1][None, :, :]
newSij = (1e21)*Sij # multiply S_ij by 1e12
Nij = sigma125[:, None, None] * id_mat[None, :, :]
newNij = (1e21)*Nij
# Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
# Sij.shape = (40, 3072, 3072)
Cij = newSij + newNij
#invCij = np.linalg.inv(Cij)
logdetC = np.linalg.slogdet(Cij) # returns sign and determinant; use logdetC[1]
# model_fit_terms = m^T C^-1 m
#
# model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) )
# for i in range(invCij.shape[0])])
#
model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ])
print "m^T c^-1 m terms are"
print model_fit_terms
print "******************"
print "logdetC terms are"
print logdetC[1]
print "******************"
print "Npix2pi term is"
print Npix2pi
""""
OUTPUT:
m^T c^-1 m terms are
[ 9.45804805e+06 1.01483017e+07 1.10268075e+07 1.20673659e+07
1.32689374e+07 1.46133094e+07 1.61825120e+07 1.79600153e+07
1.99559126e+07 2.22256705e+07 2.47646632e+07 2.76424156e+07
3.08855750e+07 3.45339903e+07 3.86587049e+07 4.33120503e+07
4.85288998e+07 5.44131447e+07 6.10398347e+07 6.84991954e+07
7.68939251e+07 8.63276615e+07 9.69909206e+07 1.08947713e+08
1.22424263e+08 1.37576408e+08 1.54634394e+08 1.73828119e+08
1.95434208e+08 2.19747753e+08 2.47115351e+08 2.77910405e+08
3.12543652e+08 3.51552581e+08 3.95426692e+08 4.44811316e+08
5.00390442e+08 5.62931595e+08 6.33312474e+08 7.12509516e+08]
******************
logdetC terms are
[ 234.83595269 -125.48871852 -486.73629436 -848.73844692
-1211.04772189 -1573.36162977 -1936.79205007 -2300.35381787
-2662.82613924 -3026.67180777 -3389.29604324 -3752.32106286
-4115.94072189 -4479.40823086 -4842.79499263 -5205.87784784
-5569.53277853 -5933.18379399 -6296.59554129 -6660.37091899
-7023.75150763 -7387.0944891 -7750.90452509 -8114.35002254
-8477.84033831 -8841.49715615 -9205.03206405 -9568.51393667
-9932.08545533 -10295.68073606 -10659.20465258 -11022.79171886
-11386.33372681 -11749.98089605 -12113.47607432 -12477.02136508
-12840.61843452 -13204.16794796 -13567.75960538 -13931.32768536]
******************
Npix2pi term is
19301.9452637
**********************************************
**********************************************
SECOND, REDUCE NOISE TO E-22 TO E-24, FIND A PEAK IN LF
CODE:
""""
vary_x_samples123 = np.logspace(-8, -12, num=40) #num = 40
sigma123 = np.logspace(-22, -24, num=40)
#
# Plot noise sigma from e-22 to e-24
#
Sij = vary_x_samples123 [:, None, None] * norm_matrix[1][None, :, :]
newSij = (1e22)*Sij # multiply S_ij by 1e12
Nij = sigma123[:, None, None] * id_mat[None, :, :]
newNij = (1e22)*Nij
# Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
# Sij.shape = (40, 3072, 3072)
Cij = newSij + newNij
#invCij = np.linalg.inv(Cij)
logdetC = np.linalg.slogdet(Cij) # returns sign and determinant; use logdetC[1]
# model_fit_terms = m^T C^-1 m
#
# model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) )
# for i in range(invCij.shape[0])])
#
model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ])
print "m^T c^-1 m terms are"
print model_fit_terms
print "******************"
print "logdetC terms are"
print logdetC[1]
print "******************"
print "Npix2pi term is"
print Npix2pi
""""
OUTPUT
m^T c^-1 m terms are
[ 1.16534465e+06 5.03576814e+06 1.18842066e+06 3.78456972e+06
2.07230699e+06 1.16453569e+06 -7.41348447e+05 -3.92254528e+06
-9.93526825e+06 -2.29984553e+07 -5.63527744e+07 -3.44097793e+08
2.13719373e+08 1.10203181e+08 9.28265870e+07 8.26463737e+07
8.06037248e+07 8.18050050e+07 8.54633954e+07 8.98256378e+07
9.79018046e+07 1.05962829e+08 1.15393283e+08 1.26566259e+08
1.39139646e+08 1.53686968e+08 1.70516270e+08 1.89068539e+08
2.10547663e+08 2.34462101e+08 2.61421501e+08 2.92204648e+08
3.26641827e+08 3.65300144e+08 4.09114198e+08 4.58308672e+08
5.13744893e+08 5.76396127e+08 6.46599075e+08 7.25782771e+08]
******************
logdetC terms are
[ -410.62898964 -672.91913446 -926.22503662 -1169.34164618
-1445.18261905 -1738.08384536 -2079.4620371 -2399.02656433
-2735.89491794 -3080.90524843 -3433.09248454 -3787.30516381
-4142.67275289 -4502.70581932 -4854.37183367 -5211.64589642
-5574.580628 -5931.36044619 -6295.41523156 -6660.63863576
-7019.91756434 -7381.27969271 -7740.57086979 -8102.90872247
-8466.51398412 -8828.15412569 -9191.80550759 -9554.54294936
-9917.99595215 -10281.00471672 -10644.32030792 -11007.88053013
-11371.45867768 -11734.62772477 -12098.20386412 -12461.47921611
-12824.82116055 -13188.76508825 -13552.04849419 -13915.45159908]
******************
Npix2pi term is
19301.9452637
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment