sugar.py 2.9 KB
 Theo Steininger committed Apr 13, 2017 1 2 3 4 5 6 7 8 9 10 11 12 ``````# 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, see . `````` Theo Steininger committed May 24, 2017 13 14 15 16 17 ``````# # Copyright(C) 2013-2017 Max-Planck-Society # # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # and financially supported by the Studienstiftung des deutschen Volkes. `````` theos committed Oct 12, 2016 18 19 20 21 `````` from nifty import PowerSpace,\ Field,\ DiagonalOperator,\ `````` Jakob Knollmueller committed May 17, 2017 22 `````` sqrt `````` theos committed Oct 12, 2016 23 24 25 26 `````` __all__ = ['create_power_operator'] `````` Jakob Knollmueller committed May 23, 2017 27 ``````def create_power_operator(domain, power_spectrum, power_domain=None, dtype=None, `````` Theo Steininger committed Mar 16, 2017 28 `````` distribution_strategy='not'): `````` Theo Steininger committed May 22, 2017 29 `````` """ Creates a diagonal operator with the given power spectrum. `````` theos committed Oct 12, 2016 30 `````` `````` Theo Steininger committed May 23, 2017 31 `````` Constructs a diagonal operator that lives over the specified domain. `````` theos committed Oct 12, 2016 32 `````` `````` Theo Steininger committed May 22, 2017 33 34 35 `````` Parameters ---------- domain : DomainObject `````` Theo Steininger committed May 23, 2017 36 `````` Domain over which the power operator shall live. `````` Theo Steininger committed May 22, 2017 37 `````` power_spectrum : (array-like, method) `````` Theo Steininger committed May 22, 2017 38 39 `````` An array-like object, or a method that implements the square root of a power spectrum as a function of k. `````` Theo Steininger committed May 22, 2017 40 `````` dtype : type *optional* `````` Theo Steininger committed May 22, 2017 41 `````` dtype that the field holding the power spectrum shall use `````` Theo Steininger committed May 22, 2017 42 43 44 `````` (default : None). if dtype == None: the dtype of `power_spectrum` will be used. distribution_strategy : string *optional* `````` Theo Steininger committed May 22, 2017 45 `````` Distributed strategy to be used by the underlying d2o objects. `````` Theo Steininger committed May 22, 2017 46 47 `````` (default : 'not') `````` Theo Steininger committed May 22, 2017 48 49 `````` Returns ------- `````` Theo Steininger committed May 22, 2017 50 `````` DiagonalOperator : An operator that implements the given power spectrum. `````` theos committed Oct 12, 2016 51 `````` `````` Theo Steininger committed May 22, 2017 52 `````` """ `````` theos committed Oct 12, 2016 53 `````` `````` theos committed Oct 15, 2016 54 55 `````` power_domain = PowerSpace(domain, distribution_strategy=distribution_strategy) `````` Martin Reinecke committed May 18, 2017 56 `````` fp = Field(power_domain, val=power_spectrum, dtype=dtype, `````` theos committed Oct 12, 2016 57 `````` distribution_strategy=distribution_strategy) `````` Theo Steininger committed May 19, 2017 58 `````` fp **= 2 `````` Pumpe, Daniel (dpumpe) committed May 11, 2017 59 `````` f = fp.power_synthesize(mean=1, std=0, real_signal=False) `````` Jakob Knollmueller committed May 17, 2017 60 `````` `````` Martin Reinecke committed May 18, 2017 61 `````` return DiagonalOperator(domain, diagonal=f, bare=True) `````` Jakob Knollmueller committed May 17, 2017 62 63 64 65 66 `````` def generate_posterior_sample(mean, covariance): S = covariance.S R = covariance.R N = covariance.N `````` Jakob Knollmueller committed May 23, 2017 67 68 69 70 `````` power = sqrt(S.diagonal().power_analyze()) mock_signal = power.power_synthesize(real_signal=True) `````` Jakob Knollmueller committed Jun 01, 2017 71 `````` noise = N.diagonal(bare=True).val `````` Jakob Knollmueller committed May 23, 2017 72 `````` `````` Jakob Knollmueller committed May 17, 2017 73 74 `````` mock_noise = Field.from_random(random_type="normal", domain=N.domain, std = sqrt(noise), dtype = noise.dtype) `````` Jakob Knollmueller committed Jun 01, 2017 75 `````` mock_data = R(mock_signal) + mock_noise `````` Jakob Knollmueller committed May 17, 2017 76 `````` `````` Jakob Knollmueller committed Jun 01, 2017 77 `````` mock_j = R.adjoint_times(N.inverse_times(mock_data)) `````` Jakob Knollmueller committed May 17, 2017 78 79 80 81 `````` mock_m = covariance.inverse_times(mock_j) sample = mock_signal - mock_m + mean return sample ``````