diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py index 7866f04d3e290f71b5457effbb2eaf69c3f7362e..75875d0eb89369a3382c9a58e5516cfb76deb354 100644 --- a/demos/critical_filtering.py +++ b/demos/critical_filtering.py @@ -7,12 +7,12 @@ np.random.seed(42) if __name__ == "__main__": # Set up position space - s_space = ift.RGSpace([128, 128]) - # s_space = ift.HPSpace(32) + #s_space = ift.RGSpace([128, 128]) + s_space = ift.HPSpace(32) # Define harmonic transformation and associated harmonic space h_space = s_space.get_default_codomain() - fft = ift.HarmonicTransformOperator(h_space, s_space) + HT = ift.HarmonicTransformOperator(h_space, s_space) # Set up power space p_space = ift.PowerSpace(h_space, @@ -34,10 +34,10 @@ if __name__ == "__main__": # Instrument._diagonal.val[64:512-64, 64:512-64] = 0 # Add a harmonic transformation to the instrument - R = Instrument*fft + R = Instrument*HT noise = 1. - N = ift.DiagonalOperator(ift.Field.full(s_space, noise).weight(1)) + N = ift.ScalingOperator(noise, s_space) n = ift.Field.from_random(domain=s_space, random_type='normal', std=np.sqrt(noise), mean=0) @@ -48,7 +48,7 @@ if __name__ == "__main__": j = R.adjoint_times(N.inverse_times(d)) realized_power = ift.log(ift.power_analyze(sh, binbounds=p_space.binbounds)) - data_power = ift.log(ift.power_analyze(fft.adjoint_times(d), + data_power = ift.log(ift.power_analyze(HT.adjoint_times(d), binbounds=p_space.binbounds)) d_data = d.val ift.plot(d, name="data.png") @@ -91,4 +91,4 @@ if __name__ == "__main__": # Plot current estimate ift.dobj.mprint(i) if i % 50 == 0: - ift.plot(fft(m0), name='map.png') + ift.plot(HT(m0), name='map.png') diff --git a/demos/log_normal_wiener_filter.py b/demos/log_normal_wiener_filter.py index a1a4ccac8346e2691ae834dec889a51a86f9630c..8f114ca73624b408ea9de07c1a0a0408a3def830 100644 --- a/demos/log_normal_wiener_filter.py +++ b/demos/log_normal_wiener_filter.py @@ -20,7 +20,7 @@ if __name__ == "__main__": # signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels) signal_space = ift.HPSpace(16) harmonic_space = signal_space.get_default_codomain() - ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) + HT = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) power_space = ift.PowerSpace(harmonic_space) # Creating the mock signal @@ -35,15 +35,14 @@ if __name__ == "__main__": # mask.val[N10*5:N10*9, N10*5:N10*9] = 0. R = ift.GeometryRemover(signal_space) R = R*ift.DiagonalOperator(mask) - R = R*ht + R = R*HT R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma) data_domain = R.target[0] # Setting up the noise covariance and drawing a random noise realization noiseless_data = R(mock_signal) noise_amplitude = noiseless_data.val.std()/signal_to_noise - N = ift.DiagonalOperator( - ift.Field.full(data_domain, noise_amplitude**2)) + N = ift.ScalingOperator(noise_amplitude**2, data_domain) noise = ift.Field.from_random( domain=data_domain, random_type='normal', std=noise_amplitude, mean=0) @@ -61,11 +60,11 @@ if __name__ == "__main__": #minimizer = ift.SteepestDescent(controller=ctrl2) me = minimizer(energy) - m = ht(me[0].position) + m = HT(me[0].position) # Plotting plotdict = {"colormap": "Planck-like"} - ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict) + ift.plot(HT(mock_signal), name="mock_signal.png", **plotdict) logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape) ift.plot(ift.Field(signal_space, val=ift.dobj.from_global_data(logdata)), name="log_of_data.png", **plotdict) @@ -75,8 +74,8 @@ if __name__ == "__main__": class Proby(ift.DiagonalProberMixin, ift.Prober): pass proby = Proby(signal_space, probe_count=1) - proby(lambda z: ht(me2[0].curvature.inverse_times(ht.adjoint_times(z)))) + proby(lambda z: HT(me2[0].curvature.inverse_times(HT.adjoint_times(z)))) sm = ift.FFTSmoothingOperator(signal_space, sigma=0.02) - variance = sm(proby.diagonal.weight(-1)) + variance = sm(proby.diagonal.weigHT(-1)) ift.plot(variance, name='variance.png', **plotdict) diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py index 501f6f0ef1f059114eca0b1c415ebf38eee30ea9..2627fdc698fff114fc42d0bf9d89914551b833e4 100644 --- a/demos/nonlinear_critical_filter.py +++ b/demos/nonlinear_critical_filter.py @@ -66,8 +66,7 @@ if __name__ == "__main__": true_sky = nonlinearity(HT(power*sh)) noiseless_data = MeasurementOperator(true_sky) noise_amplitude = noiseless_data.val.std()*noise_level - N = ift.DiagonalOperator( - ift.Field.full(d_space, noise_amplitude**2)) + N = ift.ScalingOperator(noise_amplitude**2, d_space) n = ift.Field.from_random( domain=d_space, random_type='normal', std=noise_amplitude, mean=0) diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py index 5c690f418a78b2e37aa3e62de486e537c842b661..faf60daf111bad7cc591341728e6c99ffeab35ee 100644 --- a/demos/paper_demos/cartesian_wiener_filter.py +++ b/demos/paper_demos/cartesian_wiener_filter.py @@ -83,7 +83,6 @@ if __name__ == "__main__": noiseless_data = R(mock_signal) noise_amplitude = noiseless_data.val.std()/signal_to_noise # Setting up the noise covariance and drawing a random noise realization - #ndiag = ift.Field.full(data_domain, noise_amplitude**2) N = ift.ScalingOperator(noise_amplitude**2, data_domain) noise = ift.Field.from_random( domain=data_domain, random_type='normal', @@ -100,8 +99,7 @@ if __name__ == "__main__": m_k = wiener_curvature.inverse_times(j) m = ht(m_k) - plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", - "colormap": "Planck-like"} + plotdict = {"colormap": "Planck-like"} plot_space = ift.RGSpace((N_pixels_1, N_pixels_2)) ift.plot(ift.Field(plot_space,val=ht(mock_signal).val), name='mock_signal.png', **plotdict) @@ -110,3 +108,4 @@ if __name__ == "__main__": # sampling the uncertainty map mean, variance = ift.probe_with_posterior_samples(wiener_curvature, m_k, ht, 10) ift.plot(ift.Field(plot_space, val=ift.sqrt(variance).val), name="uncertainty.png", **plotdict) + ift.plot(ift.Field(plot_space, val=mean.val), name="posterior_mean.png", **plotdict) diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py index 2d48e6c9b11d99b83436c19e62b6c9cec4f47f4c..d36218b185ac68e19c972aaebaa3ddcf56e48fc5 100644 --- a/demos/paper_demos/wiener_filter.py +++ b/demos/paper_demos/wiener_filter.py @@ -44,8 +44,7 @@ if __name__ == "__main__": noiseless_data = R(mock_signal) noise_amplitude = noiseless_data.val.std()/signal_to_noise # Setting up the noise covariance and drawing a random noise realization - ndiag = ift.Field.full(data_domain, noise_amplitude**2) - N = ift.DiagonalOperator(ndiag) + N = ift.ScalingOperator(noise_amplitude**2, data_domain) noise = ift.Field.from_random( domain=data_domain, random_type='normal', std=noise_amplitude, mean=0) @@ -60,20 +59,13 @@ if __name__ == "__main__": m_k = wiener_curvature.inverse_times(j) m = ht(m_k) - plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index", - "colormap": "Planck-like"} + plotdict = {"colormap": "Planck-like"} ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict) ift.plot(ift.Field(signal_space, val=data.val), name="data.png", **plotdict) ift.plot(m, name="map.png", **plotdict) # sampling the uncertainty map - sample_variance = ift.Field.zeros(signal_space) - sample_mean = ift.Field.zeros(signal_space) - n_samples = 10 - for i in range(n_samples): - sample = ht(wiener_curvature.draw_sample()) + m - sample_variance += sample**2 - sample_mean += sample - variance = sample_variance/n_samples - (sample_mean/n_samples)**2 + mean, variance = ift.probe_with_posterior_samples(wiener_curvature, m_k, ht, 5) ift.plot(ift.sqrt(variance), name="uncertainty.png", **plotdict) + ift.plot(mean, name="posterior_mean.png", **plotdict) diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py index 0aa9394aff24fb72a3911a0503bedc58e5370121..8e2a9d671b1f6cef3f295de426a629e6a1865250 100644 --- a/demos/wiener_filter_easy.py +++ b/demos/wiener_filter_easy.py @@ -43,7 +43,7 @@ if __name__ == "__main__": noiseless_data = R(sh) signal_to_noise = 1. noise_amplitude = noiseless_data.val.std()/signal_to_noise - N = ift.DiagonalOperator(ift.Field.full(s_space, noise_amplitude**2)) + N = ift.ScalingOperator(noise_amplitude**2, s_space) n = ift.Field.from_random(domain=s_space, random_type='normal', std=noise_amplitude, diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py index 0182f87cdb06eaf0ec4a1a936a573864b0633c86..a8493e3caa31476191893863240cd70bf95e184b 100644 --- a/demos/wiener_filter_via_curvature.py +++ b/demos/wiener_filter_via_curvature.py @@ -59,8 +59,7 @@ if __name__ == "__main__": noiseless_data = R(mock_signal) noise_amplitude = noiseless_data.val.std()/signal_to_noise - N = ift.DiagonalOperator( - ift.Field.full(data_domain, noise_amplitude**2)) + N = ift.ScalingOperator(noise_amplitude**2, data_domain) noise = ift.Field.from_random( domain=data_domain, random_type='normal', std=noise_amplitude, mean=0) diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py index 283400e878ef3a0b07a1e573b4ecb46e7a89e09d..ce1a2f05a6fe1d0c8ba6aa66fc9a79413358c37c 100644 --- a/demos/wiener_filter_via_hamiltonian.py +++ b/demos/wiener_filter_via_hamiltonian.py @@ -42,7 +42,7 @@ if __name__ == "__main__": noiseless_data = R(sh) signal_to_noise = 1. noise_amplitude = noiseless_data.val.std()/signal_to_noise - N = ift.DiagonalOperator(ift.Field.full(s_space, noise_amplitude**2)) + N = ift.ScalingOperator(noise_amplitude**2, s_space) n = ift.Field.from_random(domain=s_space, random_type='normal', std=noise_amplitude, diff --git a/nifty4/field.py b/nifty4/field.py index dd34401d2eb316ccf9446228777d6b4738786866..33435ee3db3ef03066b4a4c4fa1a4952d90617e5 100644 --- a/nifty4/field.py +++ b/nifty4/field.py @@ -529,13 +529,9 @@ class Field(object): return "<nifty4.Field>" def __str__(self): - minmax = [self.min(), self.max()] - mean = self.mean() return "nifty4.Field instance\n- domain = " + \ self._domain.__str__() + \ - "\n- val = " + repr(self.val) + \ - "\n - min.,max. = " + str(minmax) + \ - "\n - mean = " + str(mean) + "\n- val = " + repr(self.val) # Arithmetic functions working on Fields