replace integral/1-norm over amplitudes with 2-norm

nifty5/library/correlated_fields.py

@@ -171,13 +171,14 @@ def MultiCorrelatedField(target, amplitudes, gs_mean, gs_var, name='xi'):
         # Normalize amplitudes and apply spreader stack
         if amplitudes[i] is not None:
             # Dimension with correlation structure
-            normed_amplitude = OperatorNormalizer(amplitudes[i])
+            normed_amplitude = OperatorNormalizer(amplitudes[i], p=2)
             a_i.append(spr @ normed_amplitude)
         else:
             # No correlation structure in this dimension
             # -> constant amplitude
-            # still got to normalize, though.
-            amplitude_integral = Field.full(spr.domain, 1.).integrate()
+            # Still got to normalize, tough:
+            # Since the field contains only ones, 2-norm = sqrt(1-norm)
+            amplitude_integral = Field.full(spr.domain, 1.).integrate().sqrt()
             fct /= amplitude_integral
 
     fct = ScalingOperator(fct, pd.domain)

nifty5/operators/operator_normalizer.py

@@ -19,25 +19,27 @@ from ..field import Field
 from .simple_linear_operators import VdotOperator, WeightApplier
 
 
-def OperatorNormalizer(op):
+def OperatorNormalizer(op, p=1):
     """Wraps `op` into a shell that ensures the fields returned by
-    (shell @ op) will have an integral of one.
+    (shell @ op) will have a p-norm of one.
 
     Implements an operator analogously to
 
-        OperatorNormalizer(op) = op / op.integrate()
+        OperatorNormalizer(op) = op / |op|_p
 
     where the integration is implemented as an operator chain.
 
     Parameters
     ----------
     op : Operator
+    p  : int
+         Which p-norm to normalize under
     """
     ones_field = Field.full(op.target, 1.)
     vdot = VdotOperator(ones_field)
     wa = WeightApplier(op.target, spaces=None, power=1)
 
     # integrate as an operator
-    op_integrated = vdot.adjoint @ vdot @ wa @ op
+    op_integrated_inv = vdot.adjoint @ (vdot @ wa @ (op)**p)**(-1/p)
 
-    return op * op_integrated**(-1)
+    return op * op_integrated_inv