-
-
Save evanbiederstedt/a763375f068b05167c26 to your computer and use it in GitHub Desktop.
| Same as https://gist.github.com/evanbiederstedt/72b0035ca4b9fce830c3 | |
| except vary C3, hold sigma^2 fixed | |
| FIRST, NO PEAK CODE, i.e. scale covariance matrix by e+21 | |
| CODE | |
| """" | |
| # | |
| # Hold sigma^2 constant, vary C3 | |
| # | |
| vary_x_samples125 = np.logspace(-8, -12, num=40) # C3 parameter, vary from e-08 to e-12 | |
| sigma125 = 5e-22 # chose this sigma^2 parameter, hold constant | |
| Sij = vary_x_samples125[:, None, None] * norm_matrix[1][None, :, :] | |
| newSij = (1e21)*Sij # multiply S_ij by 1e12 | |
| Nij = sigma125 * 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 | |
| [ 17424944.59432721 25431124.82535299 20988175.41728179 18788639.9688897 | |
| 17380549.27187007 16487860.85815524 15909766.35792065 | |
| 15497517.16395973 15189182.91137147 14972008.39197645 | |
| 14793336.73669815 14664017.69568876 14568322.70447167 14492458.7918589 | |
| 14433833.23682215 14389138.7429928 14352751.6784684 14325011.10904521 | |
| 14302605.58943325 14285747.28930357 14272097.84990945 | |
| 14261432.05420008 14253089.98221404 14246570.26856775 | |
| 14241465.54843466 14237157.87713204 14234018.59879179 | |
| 14231438.17048982 14229450.48704928 14227802.58328702 | |
| 14226586.63568589 14225606.63455151 14224772.51924103 | |
| 14224178.22232257 14223689.86871848 14223295.19405429 | |
| 14222996.47245571 14222746.69923505 14222557.92339739 | |
| 14222409.72199761] | |
| ****************** | |
| logdetC terms are | |
| [-1910.79269502 -1899.49346143 -1894.31400072 -1893.9865907 -1890.0321336 | |
| -1889.55785882 -1889.79498435 -1890.52858638 -1891.04930553 -1892.56437866 | |
| -1893.64297778 -1895.06109716 -1896.61902812 -1898.36601353 -1899.78459156 | |
| -1901.39518245 -1903.05163531 -1904.67733691 -1906.32176722 -1907.95386497 | |
| -1909.59522293 -1911.21625436 -1912.86998408 -1914.50988852 -1916.15690315 | |
| -1917.81120861 -1919.462296 -1921.10832586 -1922.76238417 -1924.41349353 | |
| -1926.06744684 -1927.71856858 -1929.37104233 -1931.02466563 -1932.67642731 | |
| -1934.32805631 -1935.98147704 -1937.63429343 -1939.28755537 -1940.94054106] | |
| ****************** | |
| Npix2pi term is | |
| 19301.9452637 | |
| ************************************** | |
| ************************************* | |
| SECOND, THERE IS A PEAK CODE, i.e. scale covariance matrix by e+22 | |
| CODE | |
| """" | |
| # | |
| # Hold sigma^2 constant, vary C3 | |
| # | |
| vary_x_samples123 = np.logspace(-8, -12, num=40) # C3 parameter, vary from e-08 to e-12 | |
| sigma123 = 5e-22 # chose this sigma^2 parameter, hold constant | |
| # | |
| # 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 * 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 | |
| [ 2931868.99027531 2561423.09736563 2092876.99077943 1876397.00384942 | |
| 1740382.76555146 1650113.62054166 1589940.751404 1549446.41704495 | |
| 1519557.92056381 1496555.37291566 1479682.09032262 1466653.76019673 | |
| 1456837.51067493 1449220.46772181 1443528.73129329 1438909.39635894 | |
| 1435211.46326468 1432468.39780895 1430277.65089106 1428567.73979521 | |
| 1427238.59379346 1426165.74999515 1425328.54426065 1424671.52510419 | |
| 1424138.61076047 1423723.85218339 1423394.01406752 1423137.00097182 | |
| 1422944.56767479 1422779.35144192 1422656.21763844 1422559.03605542 | |
| 1422480.3001671 1422414.48572627 1422367.35435307 1422328.77614777 | |
| 1422298.52139999 1422275.64094104 1422255.88586292 1422241.38234343] | |
| ****************** | |
| logdetC terms are | |
| [ 5162.11916985 5174.89668376 5179.54501215 5178.57071109 5182.81051551 | |
| 5184.46873299 5183.19398907 5183.19703352 5182.37013867 5181.01067721 | |
| 5179.97913256 5178.52226416 5176.80653401 5175.16623507 5173.76193185 | |
| 5172.11691893 5170.45643259 5168.87473124 5167.23262542 5165.6055479 | |
| 5163.95696043 5162.3116691 5160.67301514 5159.02720484 5157.37926926 | |
| 5155.73010752 5154.07978174 5152.43292945 5150.77874406 5149.12856074 | |
| 5147.47438028 5145.82296635 5144.17003878 5142.51772009 5140.86502619 | |
| 5139.21290009 5137.56004832 5135.9069266 5134.25393896 5132.60101615] | |
| ****************** | |
| Npix2pi term is | |
| 19301.9452637 |
FIRST, FIND NO PEAK IN LF, NOISE PARAMETERS E-22 TO E-24
CODE, C3 = 5e-10:
""""
Hold C3 fixed, then do test again with 10% less C3 value
Use C3=5e-10
Noise is now varied from e-22 to e-24
vary_x_samples123 = 5e-10
sigma123 = np.logspace(-22, -24, num=40)
Sij = vary_x_samples123 * 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
[ 7.99112774e+06 9.15862622e+06 1.05174464e+07 1.21218540e+07
1.40481199e+07 1.63849383e+07 1.94345938e+07 2.28107290e+07
2.75426116e+07 3.39869501e+07 4.30589812e+07 5.86874289e+07
8.17059720e+07 1.64401397e+08 -3.15809360e+10 -1.55530576e+08
-7.00025652e+07 -3.69097741e+07 -2.01136741e+07 5.33705542e+06
4.54644402e+06 1.50969154e+07 3.64295994e+07 7.59843057e+07
-5.50891593e+07 1.77650620e+07 5.68860609e+07 1.24610626e+08
-1.18748462e+08 1.00196018e+08 3.26506311e+07 -5.96827008e+06
-9.43865379e+07 -2.91770949e+07 1.05304583e+08 -4.74799471e+07
-1.85701131e+07 1.86812899e+07 -4.27639327e+07 -5.13285598e+07]
logdetC terms are
[ 239.71699157 -122.71079135 -486.33852385 -849.25504606
-1212.4921126 -1576.26351061 -1944.85327619 -2305.753459
-2669.13822474 -3030.65420367 -3396.41771088 -3763.82736171
-4131.60575246 -4506.35071982 -4889.11466545 -5248.4037388
-5612.4490359 -5995.67825433 -6385.48032885 -6779.9332622
-7172.39370076 -7581.395549 -8041.60387949 -8517.11067861
-9002.83276846 -9446.49442793 -9839.83647962 -10172.92306214
-10420.76829252 -10623.50283171 -10792.24792309 -10928.10844176
-11027.19559602 -11100.06733525 -11170.65063404 -11215.66540551
-11261.93661065 -11284.43614935 -11306.60681018 -11326.08454794]
Npix2pi term is
19301.9452637
SECOND CODE,
C3 = 4.5e-10
Hold C3 fixed, then do test again with 10% less C3 value
Now use C3=4.5e-10
Noise is now varied from e-22 to e-24
vary_x_samples123 = 4.5e-10
sigma123 = np.logspace(-22, -24, num=40)
Sij = vary_x_samples123 * 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
[ 7.88230937e+06 9.00401772e+06 1.03091948e+07 1.18496346e+07
1.36733594e+07 1.58717504e+07 1.85223083e+07 2.29277071e+07
2.58774395e+07 3.11601583e+07 3.85718641e+07 4.95209561e+07
6.71862827e+07 9.81421763e+07 2.00515338e+08 1.79149266e+08
-1.57174217e+08 -7.01207170e+07 -3.93179729e+07 -2.01955132e+07
-6.98201421e+06 5.88648498e+06 1.83612716e+07 3.94038299e+07
5.39473080e+07 3.01503865e+07 -3.34339516e+07 -1.31643012e+08
-1.18775863e+08 -8.63387028e+06 2.02255477e+08 -5.89407096e+07
-1.31132724e+08 2.82823527e+07 8.34855351e+06 1.41266857e+08
1.68945314e+08 -8.34350175e+07 -4.03057940e+07 1.36690513e+09]
logdetC terms are
[ 238.87906146 -123.59978549 -486.97635053 -849.54919459
-1212.63010956 -1576.01604717 -1942.33073337 -2309.73321034
-2670.91431763 -3031.09627919 -3396.32588831 -3762.14204865
-4128.07960518 -4502.34341668 -4872.33980536 -5255.02597355
-5607.79232715 -5986.94061374 -6370.62068917 -6753.73665052
-7147.26442618 -7537.878678 -7981.39859695 -8433.45165572
-8907.83987127 -9384.80250567 -9848.92410084 -10252.021386
-10558.32732388 -10804.17249772 -11029.95815536 -11179.59272393
-11324.19494314 -11416.3809462 -11503.90735819 -11568.77486146
-11622.48104156 -11656.23132841 -11687.85082231 -11713.95931238]
Npix2pi term is
19301.9452637
We can do the same tests with holding C3 and varying the sigma^2 noise parameter:
First use C3 = 5e-10, and then try 10% less, i.e. C3=4.5e-10
FIRST CODE, fix C3 = 5e-10
""""
Hold C3 fixed, then do test again with 10% less C3 value
We use C3 = 5e-10
vary_x_samples125 = 5e-10
sigma125 = np.logspace(-21, -23, num=40)
Sij = vary_x_samples125 * 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
[ 7.18293292e+06 8.09386470e+06 9.12113596e+06 1.02809343e+07
1.15907061e+07 1.30724929e+07 1.47468537e+07 1.66377572e+07
1.87792946e+07 2.12083991e+07 2.39623566e+07 2.70868218e+07
3.06414381e+07 3.46776845e+07 3.92815846e+07 4.45663590e+07
5.05957546e+07 5.75247398e+07 6.54552543e+07 7.47090796e+07
8.54799215e+07 9.81015899e+07 1.12663772e+08 1.30393234e+08
1.51636111e+08 1.77587337e+08 2.08278721e+08 2.50099439e+08
3.03916267e+08 3.79636875e+08 4.98623636e+08 7.21005423e+08
1.10892079e+09 3.23917412e+09 -3.41985605e+09 -9.80515049e+08
-5.10193129e+08 -2.76047856e+08 -1.17087469e+08 -1.90858052e+07]
logdetC terms are
[ 226.9052794 -135.07131727 -496.90121724 -858.84168741
-1220.95683913 -1583.00135661 -1944.915235 -2306.69141473
-2668.93133159 -3030.70976015 -3393.04701702 -3754.86449285
-4117.08237295 -4478.81252165 -4840.70266103 -5203.60698497
-5565.30740991 -5927.80221985 -6289.9127427 -6652.62034049
-7016.06626191 -7378.97770017 -7739.66767859 -8103.60482614
-8468.13642146 -8832.06627575 -9201.75852318 -9558.79338163
-9920.82765381 -10286.47587538 -10655.88535624 -11026.26666575
-11395.66670283 -11762.37585804 -12132.20128441 -12499.66763487
-12883.12465843 -13266.56868972 -13654.89410426 -14040.44010189]
Npix2pi term is
19301.9452637
SECOND CODE, C3 fixed to 4.5e-10
"""""
Hold C3 fixed, then do test again with 10% less C3 value
Now use C3 = 4.5e-10
vary_x_samples125 = 5e-10
sigma125 = np.logspace(-21, -23, num=40)
Sij = vary_x_samples125 * 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
[ 7.18293292e+06 8.09386470e+06 9.12113596e+06 1.02809343e+07
1.15907061e+07 1.30724929e+07 1.47468537e+07 1.66377572e+07
1.87792946e+07 2.12083991e+07 2.39623566e+07 2.70868218e+07
3.06414381e+07 3.46776845e+07 3.92815846e+07 4.45663590e+07
5.05957546e+07 5.75247398e+07 6.54552543e+07 7.47090796e+07
8.54799215e+07 9.81015899e+07 1.12663772e+08 1.30393234e+08
1.51636111e+08 1.77587337e+08 2.08278721e+08 2.50099439e+08
3.03916267e+08 3.79636875e+08 4.98623636e+08 7.21005423e+08
1.10892079e+09 3.23917412e+09 -3.41985605e+09 -9.80515049e+08
-5.10193129e+08 -2.76047856e+08 -1.17087469e+08 -1.90858052e+07]
logdetC terms are
[ 226.9052794 -135.07131727 -496.90121724 -858.84168741
-1220.95683913 -1583.00135661 -1944.915235 -2306.69141473
-2668.93133159 -3030.70976015 -3393.04701702 -3754.86449285
-4117.08237295 -4478.81252165 -4840.70266103 -5203.60698497
-5565.30740991 -5927.80221985 -6289.9127427 -6652.62034049
-7016.06626191 -7378.97770017 -7739.66767859 -8103.60482614
-8468.13642146 -8832.06627575 -9201.75852318 -9558.79338163
-9920.82765381 -10286.47587538 -10655.88535624 -11026.26666575
-11395.66670283 -11762.37585804 -12132.20128441 -12499.66763487
-12883.12465843 -13266.56868972 -13654.89410426 -14040.44010189]
Npix2pi term is
19301.9452637