sugar.py 3.79 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 `````` Jakob Knollmueller committed Jul 05, 2017 23 ``````from nifty.minimization.conjugate_gradient import ConjugateGradient `````` theos committed Oct 12, 2016 24 25 26 ``````__all__ = ['create_power_operator'] `````` Jakob Knollmueller committed Jun 06, 2017 27 ``````def create_power_operator(domain, power_spectrum, 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 Jun 03, 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 `````` `````` Jakob Knollmueller committed Jun 06, 2017 54 55 `````` if isinstance(power_spectrum, Field): power_domain = power_spectrum.domain `````` Theo Steininger committed Jul 10, 2017 56 `````` else: `````` Jakob Knollmueller committed Jun 06, 2017 57 `````` power_domain = PowerSpace(domain, `````` Theo Steininger committed Jul 10, 2017 58 `````` distribution_strategy=distribution_strategy) `````` Jakob Knollmueller committed Jun 06, 2017 59 `````` `````` Martin Reinecke committed May 18, 2017 60 `````` fp = Field(power_domain, val=power_spectrum, dtype=dtype, `````` 61 62 63 `````` distribution_strategy='not') f = fp.power_synthesize(mean=1, std=0, real_signal=False, distribution_strategy=distribution_strategy) `````` Jakob Knollmueller committed Jun 06, 2017 64 `````` f **= 2 `````` Martin Reinecke committed May 18, 2017 65 `````` return DiagonalOperator(domain, diagonal=f, bare=True) `````` Jakob Knollmueller committed May 17, 2017 66 `````` `````` Theo Steininger committed Jul 10, 2017 67 `````` `````` Theo Steininger committed Jul 12, 2017 68 69 70 ``````def generate_posterior_sample(mean, covariance): """ Generates a posterior sample from a Gaussian distribution with given mean and covariance `````` Jakob Knollmueller committed Jul 05, 2017 71 `````` `````` Theo Steininger committed Jul 12, 2017 72 73 74 `````` This method generates samples by setting up the observation and reconstruction of a mock signal in order to obtain residuals of the right correlation which are added to the given mean. `````` Jakob Knollmueller committed Jul 05, 2017 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 `````` Parameters ---------- mean : Field the mean of the posterior Gaussian distribution covariance : WienerFilterCurvature The posterior correlation structure consisting of a response operator, noise covariance and prior signal covariance Returns ------- sample : Field Returns the a sample from the Gaussian of given mean and covariance. """ `````` Jakob Knollmueller committed May 17, 2017 91 92 93 `````` S = covariance.S R = covariance.R N = covariance.N `````` Jakob Knollmueller committed Jul 05, 2017 94 `````` `````` Jakob Knollmueller committed Jun 06, 2017 95 `````` power = S.diagonal().power_analyze()**.5 `````` Jakob Knollmueller committed May 23, 2017 96 97 `````` mock_signal = power.power_synthesize(real_signal=True) `````` Jakob Knollmueller committed Jun 01, 2017 98 `````` noise = N.diagonal(bare=True).val `````` Jakob Knollmueller committed May 23, 2017 99 `````` `````` Jakob Knollmueller committed May 17, 2017 100 `````` mock_noise = Field.from_random(random_type="normal", domain=N.domain, `````` Theo Steininger committed Jul 10, 2017 101 `````` std=sqrt(noise), dtype=noise.dtype) `````` Jakob Knollmueller committed Jun 01, 2017 102 `````` mock_data = R(mock_signal) + mock_noise `````` Jakob Knollmueller committed May 17, 2017 103 `````` `````` Jakob Knollmueller committed Jun 01, 2017 104 `````` mock_j = R.adjoint_times(N.inverse_times(mock_data)) `````` Jakob Knollmueller committed May 17, 2017 105 106 107 `````` mock_m = covariance.inverse_times(mock_j) sample = mock_signal - mock_m + mean return sample``````