02 conc

mlmc/estimator.py

@@ -363,7 +363,7 @@ def estimate_domain(quantity, sample_storage, quantile=None):
     return np.min(ranges[:, 0]), np.max(ranges[:, 1])
 
 
-def estimate_n_samples_for_target_variance(target_variance, prescribe_vars, n_ops, n_levels):
+def estimate_n_samples_for_target_variance(target_variance, prescribe_vars, n_ops, n_levels, vars_l2_norm=False):
     """
     Estimate optimal number of samples for individual levels that should provide a target variance of
     resulting moment estimate.
@@ -372,19 +372,26 @@ def estimate_n_samples_for_target_variance(target_variance, prescribe_vars, n_op
     :param prescribe_vars: vars[ L, M] for all levels L and moments_fn M safe the (zeroth) constant moment with zero variance.
     :param n_ops: number of operations at each level
     :param n_levels: number of levels
+    :param vars_l2_norm: if True, then use L2 norm of moments variances for each level variance
     :return: np.array with number of optimal samples for individual levels and moments_fn, array (LxR)
     """
     vars = prescribe_vars
+    if vars_l2_norm:
+        vars = np.linalg.norm(vars, axis=1)
+
     sqrt_var_n = np.sqrt(vars.T * n_ops)  # moments_fn in rows, levels in cols
-    total = np.sum(sqrt_var_n, axis=1)  # sum over levels
+    total = np.sum(sqrt_var_n, axis=-1)  # sum over levels
+
     n_samples_estimate = np.round((sqrt_var_n / n_ops).T * total / target_variance).astype(int)  # moments_fn in cols
+
     # Limit maximal number of samples per level
     n_samples_estimate_safe = np.maximum(
         np.minimum(n_samples_estimate, vars * n_levels / target_variance), 2)
 
+    if vars_l2_norm:
+        return n_samples_estimate_safe.astype(int)
     return np.max(n_samples_estimate_safe, axis=1).astype(int)
 
-
 def calc_level_params(step_range, n_levels):
     assert step_range[0] > step_range[1]
     level_parameters = []

mlmc/random/correlated_field.py

@@ -62,7 +62,7 @@ def __init__(self, name, field=None, param_fields=[], regions=[]):
         if type(field) in [float, int]:
             self.const = field
             assert len(param_fields) == 0
-        elif isinstance(field, RandomFieldBase):
+        elif isinstance(field, RandomFieldBase) or isinstance(field, gstools.covmodel.models.CovModel):
             self.correlated_field = field
             assert len(param_fields) == 0
         else:
@@ -372,7 +372,6 @@ def _initialize(self, **kwargs):
         """
         Called after initialization in common constructor.
         """
-
         ### Attributes computed in precalculation.
         self.cov_mat = None
         # Covariance matrix (dense).
@@ -500,14 +499,14 @@ def _sample(self):
 
 class GSToolsSpatialCorrelatedField(RandomFieldBase):
 
-    def __init__(self, model, mode_no=1000, log=False, sigma=1):
+    def __init__(self, model, mode_no=1000, log=False, sigma=1, mean=0):
         """
         :param model: instance of covariance model class, which parent is gstools.covmodel.CovModel
         :param mode_no: number of Fourier modes, default: 1000 as in gstools package
         """
         self.model = model
         self.mode_no = mode_no
-        self.srf = gstools.SRF(model, mode_no=mode_no)
+        self.srf = gstools.SRF(model, mean=mean, mode_no=mode_no)
         self.mu = self.srf.mean
         self.sigma = sigma
         self.dim = model.dim