sugar.py 3.97 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 28 `````` def create_power_operator(domain, power_spectrum, dtype=None, `````` Theo Steininger committed Mar 16, 2017 29 `````` distribution_strategy='not'): `````` Theo Steininger committed May 22, 2017 30 `````` """ Creates a diagonal operator with the given power spectrum. `````` theos committed Oct 12, 2016 31 `````` `````` Theo Steininger committed May 23, 2017 32 `````` Constructs a diagonal operator that lives over the specified domain. `````` theos committed Oct 12, 2016 33 `````` `````` Theo Steininger committed May 22, 2017 34 35 36 `````` Parameters ---------- domain : DomainObject `````` Theo Steininger committed Jun 03, 2017 37 `````` Domain over which the power operator shall live. `````` Theo Steininger committed May 22, 2017 38 `````` power_spectrum : (array-like, method) `````` Theo Steininger committed May 22, 2017 39 40 `````` 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 41 `````` dtype : type *optional* `````` Theo Steininger committed May 22, 2017 42 `````` dtype that the field holding the power spectrum shall use `````` Theo Steininger committed May 22, 2017 43 44 45 `````` (default : None). if dtype == None: the dtype of `power_spectrum` will be used. distribution_strategy : string *optional* `````` Theo Steininger committed May 22, 2017 46 `````` Distributed strategy to be used by the underlying d2o objects. `````` Theo Steininger committed May 22, 2017 47 48 `````` (default : 'not') `````` Theo Steininger committed May 22, 2017 49 50 `````` Returns ------- `````` Theo Steininger committed May 22, 2017 51 `````` DiagonalOperator : An operator that implements the given power spectrum. `````` theos committed Oct 12, 2016 52 `````` `````` Theo Steininger committed May 22, 2017 53 `````` """ `````` theos committed Oct 12, 2016 54 `````` `````` Jakob Knollmueller committed Jun 06, 2017 55 56 57 58 `````` if isinstance(power_spectrum, Field): power_domain = power_spectrum.domain else : power_domain = PowerSpace(domain, `````` theos committed Oct 15, 2016 59 `````` distribution_strategy=distribution_strategy) `````` Jakob Knollmueller committed Jun 06, 2017 60 61 62 `````` `````` Martin Reinecke committed May 18, 2017 63 `````` fp = Field(power_domain, val=power_spectrum, dtype=dtype, `````` theos committed Oct 12, 2016 64 `````` distribution_strategy=distribution_strategy) `````` Pumpe, Daniel (dpumpe) committed May 11, 2017 65 `````` f = fp.power_synthesize(mean=1, std=0, real_signal=False) `````` Jakob Knollmueller committed Jun 06, 2017 66 `````` f **= 2 `````` Martin Reinecke committed May 18, 2017 67 `````` return DiagonalOperator(domain, diagonal=f, bare=True) `````` Jakob Knollmueller committed May 17, 2017 68 `````` `````` Jakob Knollmueller committed Jul 05, 2017 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 ``````def generate_posterior_sample(mean, covariance, inverter = None): """ Generates a posterior sample from a Gaussian distribution with given mean and covariance 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. 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 inverter : ConjugateGradient *optional* the conjugate gradient used to invert the curvature for the Wiener filter. default : None Returns ------- sample : Field Returns the a sample from the Gaussian of given mean and covariance. """ `````` Jakob Knollmueller committed May 17, 2017 93 94 95 `````` S = covariance.S R = covariance.R N = covariance.N `````` Jakob Knollmueller committed Jul 05, 2017 96 97 98 `````` if inverter is None: inverter = ConjugateGradient(preconditioner=S) `````` Jakob Knollmueller committed Jun 06, 2017 99 `````` power = S.diagonal().power_analyze()**.5 `````` Jakob Knollmueller committed May 23, 2017 100 101 102 `````` mock_signal = power.power_synthesize(real_signal=True) `````` Jakob Knollmueller committed Jun 01, 2017 103 `````` noise = N.diagonal(bare=True).val `````` Jakob Knollmueller committed May 23, 2017 104 `````` `````` Jakob Knollmueller committed May 17, 2017 105 106 `````` mock_noise = Field.from_random(random_type="normal", domain=N.domain, std = sqrt(noise), dtype = noise.dtype) `````` Jakob Knollmueller committed Jun 01, 2017 107 `````` mock_data = R(mock_signal) + mock_noise `````` Jakob Knollmueller committed May 17, 2017 108 `````` `````` Jakob Knollmueller committed Jun 01, 2017 109 `````` mock_j = R.adjoint_times(N.inverse_times(mock_data)) `````` Jakob Knollmueller committed May 17, 2017 110 111 112 113 `````` mock_m = covariance.inverse_times(mock_j) sample = mock_signal - mock_m + mean return sample ``````