From 5e8717005241a1e5b7e550f6517acf7a835eafca Mon Sep 17 00:00:00 2001 From: pfrank Date: Thu, 14 Nov 2019 01:02:18 +0100 Subject: [PATCH] corrections for mistake in model --- demos/find_amplitude_parameters.py | 6 +- demos/getting_started_3.py | 6 +- demos/getting_started_mf.py | 10 ++-- demos/multi_amplitudes_consistency.py | 26 +++++---- demos/newamplitudes.py | 4 +- nifty5/library/correlated_fields.py | 56 ++++++++++--------- test/test_operators/test_correlated_fields.py | 8 +-- 7 files changed, 61 insertions(+), 55 deletions(-) diff --git a/demos/find_amplitude_parameters.py b/demos/find_amplitude_parameters.py index 84a68944..d470e20b 100644 --- a/demos/find_amplitude_parameters.py +++ b/demos/find_amplitude_parameters.py @@ -28,16 +28,16 @@ def _default_pspace(dom): if __name__ == '__main__': np.random.seed(42) - fa = ift.CorrelatedFieldMaker() + fa = ift.CorrelatedFieldMaker.make(10, 0.1, '') n_samps = 20 slope_means = [-2, -3] fa.add_fluctuations(ift.RGSpace(128, 0.1), 10, 2, 1, 1e-6, 2, 1e-6, slope_means[0], 0.2, 'spatial') # fa.add_fluctuations(_default_pspace(ift.RGSpace((128, 64))), 10, 2, 1, # 1e-6, 2, 1e-6, slope_means[0], 0.2, 'spatial') - fa.add_fluctuations(ift.RGSpace(32), 10, 5, 1, 1e-6, 2, + fa.add_fluctuations(ift.RGSpace(32), 3, 5, 1, 1e-6, 2, 1e-6, slope_means[1], 1, 'freq') - correlated_field = fa.finalize(10, 0.1, '') + correlated_field = fa.finalize() amplitudes = fa.amplitudes plt.style.use('seaborn-notebook') diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py index 528fc19a..3d1e66f1 100644 --- a/demos/getting_started_3.py +++ b/demos/getting_started_3.py @@ -57,10 +57,10 @@ if __name__ == '__main__': position_space = ift.RGSpace([128, 128]) - cfmaker = ift.CorrelatedFieldMaker() + cfmaker = ift.CorrelatedFieldMaker.make(1e-3, 1e-6, '') cfmaker.add_fluctuations(position_space, - 1, 1e-2, 1, .5, .1, .5, -3, 0.5, '') - correlated_field = cfmaker.finalize(1e-3, 1e-6, '') + 1., 1e-2, 1, .5, .1, .5, -3, 0.5, '') + correlated_field = cfmaker.finalize() A = cfmaker.amplitudes[0] # Apply a nonlinearity diff --git a/demos/getting_started_mf.py b/demos/getting_started_mf.py index d4987228..3d22942f 100644 --- a/demos/getting_started_mf.py +++ b/demos/getting_started_mf.py @@ -56,7 +56,7 @@ def radial_los(n_los): if __name__ == '__main__': - np.random.seed(45) + np.random.seed(43) # Choose between random line-of-sight response (mode=0) and radial lines # of sight (mode=1) @@ -71,12 +71,12 @@ if __name__ == '__main__': sp1 = ift.RGSpace(npix1) sp2 = ift.RGSpace(npix2) - cfmaker = ift.CorrelatedFieldMaker() + cfmaker = ift.CorrelatedFieldMaker.make(1e-2, 1e-6, '') amp1 = 0.5 - cfmaker.add_fluctuations(sp1, amp1, 1e-2, 1, .1, .01, .5, -2, 1., 'amp1') - cfmaker.add_fluctuations(sp2, np.sqrt(1. - amp1**2), 1e-2, 1, .1, .01, .5, + cfmaker.add_fluctuations(sp1, 0.1, 1e-2, 1, .1, .01, .5, -2, 1., 'amp1') + cfmaker.add_fluctuations(sp2, 0.1, 1e-2, 1, .1, .01, .5, -1.5, .5, 'amp2') - correlated_field = cfmaker.finalize(1e-3, 1e-6, '') + correlated_field = cfmaker.finalize() A1 = cfmaker.amplitudes[0] A2 = cfmaker.amplitudes[1] diff --git a/demos/multi_amplitudes_consistency.py b/demos/multi_amplitudes_consistency.py index e7b44508..db128551 100644 --- a/demos/multi_amplitudes_consistency.py +++ b/demos/multi_amplitudes_consistency.py @@ -11,7 +11,7 @@ def testAmplitudesConsistency(seed, sspace): sc.add(op(s.extract(op.domain))) return sc.mean.to_global_data(), sc.var.sqrt().to_global_data() np.random.seed(seed) - offset_std = 30 + offset_std = .1 intergated_fluct_std0 = .003 intergated_fluct_std1 = 0.1 @@ -20,12 +20,12 @@ def testAmplitudesConsistency(seed, sspace): fsspace = ift.RGSpace((12,), (0.4,)) - fa = ift.CorrelatedFieldMaker() + fa = ift.CorrelatedFieldMaker.make(offset_std, 1E-8, '') fa.add_fluctuations(sspace, intergated_fluct_std0, 1E-8, 1.1, 2., 2.1, .5, -2, 1., 'spatial') fa.add_fluctuations(fsspace, intergated_fluct_std1, 1E-8, 3.1, 1., .5, .1, -4, 1., 'freq') - op = fa.finalize(offset_std, 1E-8, '') + op = fa.finalize() samples = [ift.from_random('normal',op.domain) for _ in range(nsam)] tot_flm, _ = stats(fa.total_fluctuation,samples) @@ -44,13 +44,6 @@ def testAmplitudesConsistency(seed, sspace): sl_fluct_space = fa.slice_fluctuation_realized(sams,0) sl_fluct_freq = fa.slice_fluctuation_realized(sams,1) - np.testing.assert_allclose(offset_std, zm_std_mean, rtol=0.5) - np.testing.assert_allclose(intergated_fluct_std0, fluct_space, rtol=0.5) - np.testing.assert_allclose(intergated_fluct_std1, fluct_freq, rtol=0.5) - np.testing.assert_allclose(tot_flm, fluct_total, rtol=0.5) - np.testing.assert_allclose(slice_fluct_std0, sl_fluct_space, rtol=0.5) - np.testing.assert_allclose(slice_fluct_std1, sl_fluct_freq, rtol=0.5) - print("Expected offset Std: " + str(offset_std)) print("Estimated offset Std: " + str(zm_std_mean)) @@ -73,14 +66,23 @@ def testAmplitudesConsistency(seed, sspace): print("Expected total fluct. Std: " + str(tot_flm)) print("Estimated total fluct. Std: " + str(fluct_total)) - fa = ift.CorrelatedFieldMaker() + + np.testing.assert_allclose(offset_std, zm_std_mean, rtol=0.5) + np.testing.assert_allclose(intergated_fluct_std0, fluct_space, rtol=0.5) + np.testing.assert_allclose(intergated_fluct_std1, fluct_freq, rtol=0.5) + np.testing.assert_allclose(tot_flm, fluct_total, rtol=0.5) + np.testing.assert_allclose(slice_fluct_std0, sl_fluct_space, rtol=0.5) + np.testing.assert_allclose(slice_fluct_std1, sl_fluct_freq, rtol=0.5) + + + fa = ift.CorrelatedFieldMaker.make(offset_std, .1, '') fa.add_fluctuations(fsspace, intergated_fluct_std1, 1., 3.1, 1., .5, .1, -4, 1., 'freq') m = 3. x = fa.moment_slice_to_average(m) fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5, -2, 1., 'spatial', 0) - op = fa.finalize(offset_std, .1, '') + op = fa.finalize() em, estd = stats(fa.slice_fluctuation(0),samples) print("Forced slice fluct. space Std: "+str(m)) print("Expected slice fluct. Std: " + str(em)) diff --git a/demos/newamplitudes.py b/demos/newamplitudes.py index 507cdbb5..b0164ab9 100644 --- a/demos/newamplitudes.py +++ b/demos/newamplitudes.py @@ -4,9 +4,9 @@ np.random.seed(42) sspace = ift.RGSpace((128,)) -fa = ift.CorrelatedFieldMaker() +fa = ift.CorrelatedFieldMaker.make(10, 0.1, '') fa.add_fluctuations(sspace, 10, 2, 1, 1e-6, 2, 1e-6, -2, 1e-6, 'spatial') -op = fa.finalize(10, 0.1, '') +op = fa.finalize() A = fa.amplitudes[0] cstpos = ift.from_random('normal', op.domain) diff --git a/nifty5/library/correlated_fields.py b/nifty5/library/correlated_fields.py index 6659f822..d4db894d 100644 --- a/nifty5/library/correlated_fields.py +++ b/nifty5/library/correlated_fields.py @@ -255,10 +255,23 @@ class _Amplitude(Operator): class CorrelatedFieldMaker: - def __init__(self): + def __init__(self, amplitude_offset, prefix): self._a = [] - self._azm = None self._position_spaces = [] + + self._azm = amplitude_offset + self._prefix = prefix + + @staticmethod + def make(offset_amplitude_mean, offset_amplitude_stddev, prefix): + offset_amplitude_stddev = float(offset_amplitude_stddev) + offset_amplitude_mean = float(offset_amplitude_mean) + assert offset_amplitude_stddev > 0 + assert offset_amplitude_mean > 0 + zm = _LognormalMomentMatching(offset_amplitude_mean, + offset_amplitude_stddev, + prefix + 'zeromode') + return CorrelatedFieldMaker(zm, prefix) def add_fluctuations(self, position_space, @@ -292,6 +305,7 @@ class CorrelatedFieldMaker: fluct = _LognormalMomentMatching(fluctuations_mean, fluctuations_stddev, prefix + 'fluctuations') + fluct = fluct*self._azm.one_over() flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev, prefix + 'flexibility') asp = _LognormalMomentMatching(asperity_mean, asperity_stddev, @@ -309,14 +323,12 @@ class CorrelatedFieldMaker: def finalize_from_op(self, zeromode, prefix=''): assert isinstance(zeromode, Operator) - self._azm = zeromode hspace = makeDomain([dd.get_default_codomain() for dd in self._position_spaces]) foo = np.ones(hspace.shape) zeroind = len(hspace.shape)*(0,) foo[zeroind] = 0 - azm = Adder(from_global_data(hspace, foo)) @ ValueInserter( - hspace, zeroind) @ zeromode + azm = VdotOperator(full(hspace,1.)).adjoint @ zeromode n_amplitudes = len(self._a) ht = HarmonicTransformOperator(hspace, self._position_spaces[0], @@ -340,25 +352,16 @@ class CorrelatedFieldMaker: return ht(azm*(pd @ a)*ducktape(hspace, None, prefix + 'xi')) def finalize(self, - offset_amplitude_mean, - offset_amplitude_stddev, - prefix='', offset=None, prior_info=100): """ offset vs zeromode: volume factor """ - offset_amplitude_stddev = float(offset_amplitude_stddev) - offset_amplitude_mean = float(offset_amplitude_mean) - assert offset_amplitude_stddev > 0 - assert offset_amplitude_mean > 0 if offset is not None: raise NotImplementedError offset = float(offset) - azm = _LognormalMomentMatching(offset_amplitude_mean, - offset_amplitude_stddev, - prefix + 'zeromode') - op = self.finalize_from_op(azm, prefix) + + op = self.finalize_from_op(self._azm, self._prefix) if prior_info > 0: from ..sugar import from_random samps = [ @@ -386,12 +389,13 @@ class CorrelatedFieldMaker: def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000): fluctuations_slice_mean = float(fluctuations_slice_mean) assert fluctuations_slice_mean > 0 + from ..sugar import from_random scm = 1. for a in self._a: - m, std = a.fluctuation_amplitude.mean, a.fluctuation_amplitude.std - mu, sig = _lognormal_moments(m, std) - flm = np.exp(mu + sig*np.random.normal(size=nsamples)) - scm *= flm**2 + 1. + op = a.fluctuation_amplitude + res= np.array([op(from_random('normal',op.domain)).to_global_data() + for _ in range(nsamples)]) + scm *= res**2 + 1. return fluctuations_slice_mean/np.mean(np.sqrt(scm)) @property @@ -408,12 +412,12 @@ class CorrelatedFieldMaker: if len(self._a) == 0: raise NotImplementedError if len(self._a) == 1: - return self._a[0].fluctuation_amplitude + return self.average_fluctuation(0) q = 1. for a in self._a: fl = a.fluctuation_amplitude q = q*(Adder(full(fl.target, 1.)) @ fl**2) - return (Adder(full(q.target, -1.)) @ q).sqrt() + return (Adder(full(q.target, -1.)) @ q).sqrt() * self._azm def slice_fluctuation(self, space): """Returns operator which acts on prior or posterior samples""" @@ -421,7 +425,7 @@ class CorrelatedFieldMaker: raise NotImplementedError assert space < len(self._a) if len(self._a) == 1: - return self._a[0].fluctuation_amplitude + return self.average_fluctuation(0) q = 1. for j in range(len(self._a)): fl = self._a[j].fluctuation_amplitude @@ -429,7 +433,7 @@ class CorrelatedFieldMaker: q = q*fl**2 else: q = q*(Adder(full(fl.target, 1.)) @ fl**2) - return q.sqrt() + return q.sqrt() * self._azm def average_fluctuation(self, space): """Returns operator which acts on prior or posterior samples""" @@ -437,8 +441,8 @@ class CorrelatedFieldMaker: raise NotImplementedError assert space < len(self._a) if len(self._a) == 1: - return self._a[0].fluctuation_amplitude - return self._a[space].fluctuation_amplitude + return self._a[0].fluctuation_amplitude*self._azm + return self._a[space].fluctuation_amplitude*self._azm @staticmethod def offset_amplitude_realized(samples): diff --git a/test/test_operators/test_correlated_fields.py b/test/test_operators/test_correlated_fields.py index 4beb6082..fc761fb6 100644 --- a/test/test_operators/test_correlated_fields.py +++ b/test/test_operators/test_correlated_fields.py @@ -42,12 +42,12 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std): fsspace = ift.RGSpace((12,), (0.4,)) - fa = ift.CorrelatedFieldMaker() + fa = ift.CorrelatedFieldMaker.make(offset_std, 1E-8, '') fa.add_fluctuations(sspace, Astds[0], 1E-8, 1.1, 2., 2.1, .5, -2, 1., 'spatial') fa.add_fluctuations(fsspace, Astds[1], 1E-8, 3.1, 1., .5, .1, -4, 1., 'freq') - op = fa.finalize(offset_std, 1E-8, '') + op = fa.finalize() samples = [ift.from_random('normal',op.domain) for _ in range(nsam)] tot_flm, _ = stats(fa.total_fluctuation,samples) @@ -73,14 +73,14 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std): assert_allclose(slice_fluct_std0, sl_fluct_space, rtol=0.5) assert_allclose(slice_fluct_std1, sl_fluct_freq, rtol=0.5) - fa = ift.CorrelatedFieldMaker() + fa = ift.CorrelatedFieldMaker.make(offset_std, .1, '') fa.add_fluctuations(fsspace, Astds[1], 1., 3.1, 1., .5, .1, -4, 1., 'freq') m = 3. x = fa.moment_slice_to_average(m) fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5, -2, 1., 'spatial', 0) - op = fa.finalize(offset_std, .1, '') + op = fa.finalize() em, estd = stats(fa.slice_fluctuation(0),samples) assert_allclose(m, em, rtol=0.5) -- GitLab