consistency tests for amplitude model hyperparameters

demos/multi_amplitudes_consistency.py

+import nifty5 as ift
+import numpy as np
+#np.random.seed(42)
+
+offset_std = 40
+intergated_fluct_std0 = 10.
+intergated_fluct_std1 = 2.
+
+sspace = ift.RGSpace((32,64),(1.1,0.3))
+hspace = sspace.get_default_codomain()
+target0 = ift.PowerSpace(hspace)
+
+fsspace = ift.RGSpace((12,),(0.4,))
+fhspace = fsspace.get_default_codomain()
+target1 = ift.PowerSpace(fhspace)
+
+fa = ift.CorrelatedFieldMaker()
+fa.add_fluctuations(target0, intergated_fluct_std0, 1E-8,
+                    1.1, 2., 2.1, .5, -2, 1., 'spatial')
+fa.add_fluctuations(target1, intergated_fluct_std1, 1E-8,
+                    3.1, 1., .5, .1, -4, 1., 'freq')
+op = fa.finalize(offset_std, 1E-8, '')
+
+flucts = [intergated_fluct_std0, intergated_fluct_std1]
+tot_flm,totflsig = fa.effective_total_fluctuation(flucts,[1E-8,1E-8])
+
+space = op.target
+totaltoalvol = 1.
+for s in space[:]:
+    totaltoalvol *= s.total_volume
+
+nsam = 1000
+
+zm_std_mean = 0.
+fluct_space = 0.
+fluct_freq = 0.
+fluct_total = 0.
+
+for i in range(nsam):
+    x = ift.from_random('normal',op.domain)
+    res = op(x)
+
+    zm = res.integrate()/totaltoalvol
+    zm2 = res.mean()
+    
+    fl = ((res-zm)**2).integrate() / totaltoalvol
+    
+    zm_std_mean += zm**2
+    fluct_total += fl
+
+    r = res.integrate(1)/fsspace.total_volume
+    r0 = r.integrate()/sspace.total_volume
+    tm = ((r-r0)**2).integrate() / sspace.total_volume
+    fluct_space += tm
+    
+    fr = res.integrate(0)/sspace.total_volume
+    fr0 = fr.integrate()/fsspace.total_volume
+    ftm = ((fr-fr0)**2).integrate() / fsspace.total_volume
+    fluct_freq += ftm
+    
+
+fluct_total = np.sqrt(fluct_total/nsam)
+fluct_space = np.sqrt(fluct_space/nsam)
+fluct_freq = np.sqrt(fluct_freq/nsam)
+zm_std_mean = np.sqrt(zm_std_mean/nsam)
+
+print("Expected  offset Std: "+str(offset_std))
+print("Estimated offset Std: "+str(zm_std_mean))
+
+print("Expected  integrated fluct. space Std: "+str(intergated_fluct_std0))
+print("Estimated integrated fluct. space Std: "+str(fluct_space))
+
+print("Expected  integrated fluct. frequency Std: "+str(intergated_fluct_std1))
+print("Estimated integrated fluct. frequency Std: "+str(fluct_freq))
+
+print("Expected  total fluct. Std: "+str(tot_flm))
+print("Estimated total fluct. Std: "+str(fluct_total))
\ No newline at end of file

nifty5/library/correlated_fields.py

@@ -24,7 +24,6 @@ from ..domains.power_space import PowerSpace
 from ..domains.unstructured_domain import UnstructuredDomain
 from ..extra import check_jacobian_consistency, consistency_check
 from ..field import Field
-from ..multi_domain import MultiDomain
 from ..operators.adder import Adder
 from ..operators.contraction_operator import ContractionOperator
 from ..operators.distributors import PowerDistributor
@@ -36,12 +35,16 @@ from ..operators.simple_linear_operators import VdotOperator, ducktape
 from ..operators.value_inserter import ValueInserter
 from ..sugar import from_global_data, from_random, full, makeDomain
 
-
-def _lognormal_moment_matching(mean, sig, key):
-    mean, sig, key = float(mean), float(sig), str(key)
+def _lognormal_moments(mean, sig):
+    mean, sig = float(mean), float(sig)
     assert sig > 0
     logsig = np.sqrt(np.log((sig/mean)**2 + 1))
     logmean = np.log(mean) - logsig**2/2
+    return logmean, logsig
+
+def _lognormal_moment_matching(mean, sig, key):
+    key = str(key)
+    logmean, logsig = _lognormal_moments(mean, sig)
     return _normal(logmean, logsig, key).exp()
 
 
@@ -304,3 +307,18 @@ class CorrelatedFieldMaker:
     @property
     def amplitudes(self):
         return self._amplitudes
+
+    def effective_total_fluctuation(self,fluctuations_means,
+                                    fluctuations_stddevs,
+                                    nsamples = 100):
+        namps = len(fluctuations_means)
+        xis = np.random.normal(size=namps*nsamples).reshape((namps,nsamples))
+        
+        q=np.ones(nsamples)
+        for i in range(len(fluctuations_means)):
+            m, sig = _lognormal_moments(fluctuations_means[i],
+                                        fluctuations_stddevs[i])
+            f = np.exp(m + sig*xis[i])
+            q *= (1.+ f**2)
+        q = np.sqrt(q-1.)
+        return np.mean(q), np.std(q)