Commit 5e871700 authored by Philipp Frank's avatar Philipp Frank

corrections for mistake in model

parent 05dde3a1
Pipeline #63622 passed with stages
in 22 minutes and 45 seconds
...@@ -28,16 +28,16 @@ def _default_pspace(dom): ...@@ -28,16 +28,16 @@ def _default_pspace(dom):
if __name__ == '__main__': if __name__ == '__main__':
np.random.seed(42) np.random.seed(42)
fa = ift.CorrelatedFieldMaker() fa = ift.CorrelatedFieldMaker.make(10, 0.1, '')
n_samps = 20 n_samps = 20
slope_means = [-2, -3] slope_means = [-2, -3]
fa.add_fluctuations(ift.RGSpace(128, 0.1), 10, 2, 1, 1e-6, fa.add_fluctuations(ift.RGSpace(128, 0.1), 10, 2, 1, 1e-6,
2, 1e-6, slope_means[0], 0.2, 'spatial') 2, 1e-6, slope_means[0], 0.2, 'spatial')
# fa.add_fluctuations(_default_pspace(ift.RGSpace((128, 64))), 10, 2, 1, # fa.add_fluctuations(_default_pspace(ift.RGSpace((128, 64))), 10, 2, 1,
# 1e-6, 2, 1e-6, slope_means[0], 0.2, 'spatial') # 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') 1e-6, slope_means[1], 1, 'freq')
correlated_field = fa.finalize(10, 0.1, '') correlated_field = fa.finalize()
amplitudes = fa.amplitudes amplitudes = fa.amplitudes
plt.style.use('seaborn-notebook') plt.style.use('seaborn-notebook')
......
...@@ -57,10 +57,10 @@ if __name__ == '__main__': ...@@ -57,10 +57,10 @@ if __name__ == '__main__':
position_space = ift.RGSpace([128, 128]) position_space = ift.RGSpace([128, 128])
cfmaker = ift.CorrelatedFieldMaker() cfmaker = ift.CorrelatedFieldMaker.make(1e-3, 1e-6, '')
cfmaker.add_fluctuations(position_space, cfmaker.add_fluctuations(position_space,
1, 1e-2, 1, .5, .1, .5, -3, 0.5, '') 1., 1e-2, 1, .5, .1, .5, -3, 0.5, '')
correlated_field = cfmaker.finalize(1e-3, 1e-6, '') correlated_field = cfmaker.finalize()
A = cfmaker.amplitudes[0] A = cfmaker.amplitudes[0]
# Apply a nonlinearity # Apply a nonlinearity
......
...@@ -56,7 +56,7 @@ def radial_los(n_los): ...@@ -56,7 +56,7 @@ def radial_los(n_los):
if __name__ == '__main__': if __name__ == '__main__':
np.random.seed(45) np.random.seed(43)
# Choose between random line-of-sight response (mode=0) and radial lines # Choose between random line-of-sight response (mode=0) and radial lines
# of sight (mode=1) # of sight (mode=1)
...@@ -71,12 +71,12 @@ if __name__ == '__main__': ...@@ -71,12 +71,12 @@ if __name__ == '__main__':
sp1 = ift.RGSpace(npix1) sp1 = ift.RGSpace(npix1)
sp2 = ift.RGSpace(npix2) sp2 = ift.RGSpace(npix2)
cfmaker = ift.CorrelatedFieldMaker() cfmaker = ift.CorrelatedFieldMaker.make(1e-2, 1e-6, '')
amp1 = 0.5 amp1 = 0.5
cfmaker.add_fluctuations(sp1, amp1, 1e-2, 1, .1, .01, .5, -2, 1., 'amp1') cfmaker.add_fluctuations(sp1, 0.1, 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(sp2, 0.1, 1e-2, 1, .1, .01, .5,
-1.5, .5, 'amp2') -1.5, .5, 'amp2')
correlated_field = cfmaker.finalize(1e-3, 1e-6, '') correlated_field = cfmaker.finalize()
A1 = cfmaker.amplitudes[0] A1 = cfmaker.amplitudes[0]
A2 = cfmaker.amplitudes[1] A2 = cfmaker.amplitudes[1]
......
...@@ -11,7 +11,7 @@ def testAmplitudesConsistency(seed, sspace): ...@@ -11,7 +11,7 @@ def testAmplitudesConsistency(seed, sspace):
sc.add(op(s.extract(op.domain))) sc.add(op(s.extract(op.domain)))
return sc.mean.to_global_data(), sc.var.sqrt().to_global_data() return sc.mean.to_global_data(), sc.var.sqrt().to_global_data()
np.random.seed(seed) np.random.seed(seed)
offset_std = 30 offset_std = .1
intergated_fluct_std0 = .003 intergated_fluct_std0 = .003
intergated_fluct_std1 = 0.1 intergated_fluct_std1 = 0.1
...@@ -20,12 +20,12 @@ def testAmplitudesConsistency(seed, sspace): ...@@ -20,12 +20,12 @@ def testAmplitudesConsistency(seed, sspace):
fsspace = ift.RGSpace((12,), (0.4,)) 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, fa.add_fluctuations(sspace, intergated_fluct_std0, 1E-8, 1.1, 2., 2.1, .5,
-2, 1., 'spatial') -2, 1., 'spatial')
fa.add_fluctuations(fsspace, intergated_fluct_std1, 1E-8, 3.1, 1., .5, .1, fa.add_fluctuations(fsspace, intergated_fluct_std1, 1E-8, 3.1, 1., .5, .1,
-4, 1., 'freq') -4, 1., 'freq')
op = fa.finalize(offset_std, 1E-8, '') op = fa.finalize()
samples = [ift.from_random('normal',op.domain) for _ in range(nsam)] samples = [ift.from_random('normal',op.domain) for _ in range(nsam)]
tot_flm, _ = stats(fa.total_fluctuation,samples) tot_flm, _ = stats(fa.total_fluctuation,samples)
...@@ -44,13 +44,6 @@ def testAmplitudesConsistency(seed, sspace): ...@@ -44,13 +44,6 @@ def testAmplitudesConsistency(seed, sspace):
sl_fluct_space = fa.slice_fluctuation_realized(sams,0) sl_fluct_space = fa.slice_fluctuation_realized(sams,0)
sl_fluct_freq = fa.slice_fluctuation_realized(sams,1) 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("Expected offset Std: " + str(offset_std))
print("Estimated offset Std: " + str(zm_std_mean)) print("Estimated offset Std: " + str(zm_std_mean))
...@@ -73,14 +66,23 @@ def testAmplitudesConsistency(seed, sspace): ...@@ -73,14 +66,23 @@ def testAmplitudesConsistency(seed, sspace):
print("Expected total fluct. Std: " + str(tot_flm)) print("Expected total fluct. Std: " + str(tot_flm))
print("Estimated total fluct. Std: " + str(fluct_total)) 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, fa.add_fluctuations(fsspace, intergated_fluct_std1, 1., 3.1, 1., .5, .1,
-4, 1., 'freq') -4, 1., 'freq')
m = 3. m = 3.
x = fa.moment_slice_to_average(m) x = fa.moment_slice_to_average(m)
fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5, fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5,
-2, 1., 'spatial', 0) -2, 1., 'spatial', 0)
op = fa.finalize(offset_std, .1, '') op = fa.finalize()
em, estd = stats(fa.slice_fluctuation(0),samples) em, estd = stats(fa.slice_fluctuation(0),samples)
print("Forced slice fluct. space Std: "+str(m)) print("Forced slice fluct. space Std: "+str(m))
print("Expected slice fluct. Std: " + str(em)) print("Expected slice fluct. Std: " + str(em))
......
...@@ -4,9 +4,9 @@ np.random.seed(42) ...@@ -4,9 +4,9 @@ np.random.seed(42)
sspace = ift.RGSpace((128,)) 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') 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] A = fa.amplitudes[0]
cstpos = ift.from_random('normal', op.domain) cstpos = ift.from_random('normal', op.domain)
......
...@@ -255,10 +255,23 @@ class _Amplitude(Operator): ...@@ -255,10 +255,23 @@ class _Amplitude(Operator):
class CorrelatedFieldMaker: class CorrelatedFieldMaker:
def __init__(self): def __init__(self, amplitude_offset, prefix):
self._a = [] self._a = []
self._azm = None
self._position_spaces = [] 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, def add_fluctuations(self,
position_space, position_space,
...@@ -292,6 +305,7 @@ class CorrelatedFieldMaker: ...@@ -292,6 +305,7 @@ class CorrelatedFieldMaker:
fluct = _LognormalMomentMatching(fluctuations_mean, fluct = _LognormalMomentMatching(fluctuations_mean,
fluctuations_stddev, fluctuations_stddev,
prefix + 'fluctuations') prefix + 'fluctuations')
fluct = fluct*self._azm.one_over()
flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev, flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev,
prefix + 'flexibility') prefix + 'flexibility')
asp = _LognormalMomentMatching(asperity_mean, asperity_stddev, asp = _LognormalMomentMatching(asperity_mean, asperity_stddev,
...@@ -309,14 +323,12 @@ class CorrelatedFieldMaker: ...@@ -309,14 +323,12 @@ class CorrelatedFieldMaker:
def finalize_from_op(self, zeromode, prefix=''): def finalize_from_op(self, zeromode, prefix=''):
assert isinstance(zeromode, Operator) assert isinstance(zeromode, Operator)
self._azm = zeromode
hspace = makeDomain([dd.get_default_codomain() hspace = makeDomain([dd.get_default_codomain()
for dd in self._position_spaces]) for dd in self._position_spaces])
foo = np.ones(hspace.shape) foo = np.ones(hspace.shape)
zeroind = len(hspace.shape)*(0,) zeroind = len(hspace.shape)*(0,)
foo[zeroind] = 0 foo[zeroind] = 0
azm = Adder(from_global_data(hspace, foo)) @ ValueInserter( azm = VdotOperator(full(hspace,1.)).adjoint @ zeromode
hspace, zeroind) @ zeromode
n_amplitudes = len(self._a) n_amplitudes = len(self._a)
ht = HarmonicTransformOperator(hspace, self._position_spaces[0], ht = HarmonicTransformOperator(hspace, self._position_spaces[0],
...@@ -340,25 +352,16 @@ class CorrelatedFieldMaker: ...@@ -340,25 +352,16 @@ class CorrelatedFieldMaker:
return ht(azm*(pd @ a)*ducktape(hspace, None, prefix + 'xi')) return ht(azm*(pd @ a)*ducktape(hspace, None, prefix + 'xi'))
def finalize(self, def finalize(self,
offset_amplitude_mean,
offset_amplitude_stddev,
prefix='',
offset=None, offset=None,
prior_info=100): prior_info=100):
""" """
offset vs zeromode: volume factor 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: if offset is not None:
raise NotImplementedError raise NotImplementedError
offset = float(offset) offset = float(offset)
azm = _LognormalMomentMatching(offset_amplitude_mean,
offset_amplitude_stddev, op = self.finalize_from_op(self._azm, self._prefix)
prefix + 'zeromode')
op = self.finalize_from_op(azm, prefix)
if prior_info > 0: if prior_info > 0:
from ..sugar import from_random from ..sugar import from_random
samps = [ samps = [
...@@ -386,12 +389,13 @@ class CorrelatedFieldMaker: ...@@ -386,12 +389,13 @@ class CorrelatedFieldMaker:
def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000): def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
fluctuations_slice_mean = float(fluctuations_slice_mean) fluctuations_slice_mean = float(fluctuations_slice_mean)
assert fluctuations_slice_mean > 0 assert fluctuations_slice_mean > 0
from ..sugar import from_random
scm = 1. scm = 1.
for a in self._a: for a in self._a:
m, std = a.fluctuation_amplitude.mean, a.fluctuation_amplitude.std op = a.fluctuation_amplitude
mu, sig = _lognormal_moments(m, std) res= np.array([op(from_random('normal',op.domain)).to_global_data()
flm = np.exp(mu + sig*np.random.normal(size=nsamples)) for _ in range(nsamples)])
scm *= flm**2 + 1. scm *= res**2 + 1.
return fluctuations_slice_mean/np.mean(np.sqrt(scm)) return fluctuations_slice_mean/np.mean(np.sqrt(scm))
@property @property
...@@ -408,12 +412,12 @@ class CorrelatedFieldMaker: ...@@ -408,12 +412,12 @@ class CorrelatedFieldMaker:
if len(self._a) == 0: if len(self._a) == 0:
raise NotImplementedError raise NotImplementedError
if len(self._a) == 1: if len(self._a) == 1:
return self._a[0].fluctuation_amplitude return self.average_fluctuation(0)
q = 1. q = 1.
for a in self._a: for a in self._a:
fl = a.fluctuation_amplitude fl = a.fluctuation_amplitude
q = q*(Adder(full(fl.target, 1.)) @ fl**2) 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): def slice_fluctuation(self, space):
"""Returns operator which acts on prior or posterior samples""" """Returns operator which acts on prior or posterior samples"""
...@@ -421,7 +425,7 @@ class CorrelatedFieldMaker: ...@@ -421,7 +425,7 @@ class CorrelatedFieldMaker:
raise NotImplementedError raise NotImplementedError
assert space < len(self._a) assert space < len(self._a)
if len(self._a) == 1: if len(self._a) == 1:
return self._a[0].fluctuation_amplitude return self.average_fluctuation(0)
q = 1. q = 1.
for j in range(len(self._a)): for j in range(len(self._a)):
fl = self._a[j].fluctuation_amplitude fl = self._a[j].fluctuation_amplitude
...@@ -429,7 +433,7 @@ class CorrelatedFieldMaker: ...@@ -429,7 +433,7 @@ class CorrelatedFieldMaker:
q = q*fl**2 q = q*fl**2
else: else:
q = q*(Adder(full(fl.target, 1.)) @ fl**2) q = q*(Adder(full(fl.target, 1.)) @ fl**2)
return q.sqrt() return q.sqrt() * self._azm
def average_fluctuation(self, space): def average_fluctuation(self, space):
"""Returns operator which acts on prior or posterior samples""" """Returns operator which acts on prior or posterior samples"""
...@@ -437,8 +441,8 @@ class CorrelatedFieldMaker: ...@@ -437,8 +441,8 @@ class CorrelatedFieldMaker:
raise NotImplementedError raise NotImplementedError
assert space < len(self._a) assert space < len(self._a)
if len(self._a) == 1: if len(self._a) == 1:
return self._a[0].fluctuation_amplitude return self._a[0].fluctuation_amplitude*self._azm
return self._a[space].fluctuation_amplitude return self._a[space].fluctuation_amplitude*self._azm
@staticmethod @staticmethod
def offset_amplitude_realized(samples): def offset_amplitude_realized(samples):
......
...@@ -42,12 +42,12 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std): ...@@ -42,12 +42,12 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std):
fsspace = ift.RGSpace((12,), (0.4,)) 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, fa.add_fluctuations(sspace, Astds[0], 1E-8, 1.1, 2., 2.1, .5,
-2, 1., 'spatial') -2, 1., 'spatial')
fa.add_fluctuations(fsspace, Astds[1], 1E-8, 3.1, 1., .5, .1, fa.add_fluctuations(fsspace, Astds[1], 1E-8, 3.1, 1., .5, .1,
-4, 1., 'freq') -4, 1., 'freq')
op = fa.finalize(offset_std, 1E-8, '') op = fa.finalize()
samples = [ift.from_random('normal',op.domain) for _ in range(nsam)] samples = [ift.from_random('normal',op.domain) for _ in range(nsam)]
tot_flm, _ = stats(fa.total_fluctuation,samples) tot_flm, _ = stats(fa.total_fluctuation,samples)
...@@ -73,14 +73,14 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std): ...@@ -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_std0, sl_fluct_space, rtol=0.5)
assert_allclose(slice_fluct_std1, sl_fluct_freq, 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, fa.add_fluctuations(fsspace, Astds[1], 1., 3.1, 1., .5, .1,
-4, 1., 'freq') -4, 1., 'freq')
m = 3. m = 3.
x = fa.moment_slice_to_average(m) x = fa.moment_slice_to_average(m)
fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5, fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5,
-2, 1., 'spatial', 0) -2, 1., 'spatial', 0)
op = fa.finalize(offset_std, .1, '') op = fa.finalize()
em, estd = stats(fa.slice_fluctuation(0),samples) em, estd = stats(fa.slice_fluctuation(0),samples)
assert_allclose(m, em, rtol=0.5) assert_allclose(m, em, rtol=0.5)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment