Changing generation of pseudodata

Created by: andreab1997

Currently, the pseudodata are generated again each time even one of them is negative (and not ASY) before going on. From a theoretical point of view this is strange, because it is equivalent to cut the gaussians distribution before zero. Moreover, from a performance point of view, this easily becomes a bottleneck (for example in the case of fit with thcovmat in which the shifts are bigger). We want to try doing something slightly different, namely taking the absolute value of each negative pseudodata. From a theoretical point of view this makes no sense as the actual method, but, at least, it requires only a single iteration in every case.

validphys2/src/validphys/pseudodata.py

@@ -198,38 +198,46 @@ def make_replica(
     full_mask=np.concatenate(check_positive_masks, axis=0)
     # The inner while True loop is for ensuring a positive definite
     # pseudodata replica
-    while True:
-        mult_shifts = []
-        # Prepare the per-dataset multiplicative shifts
-        for mult_uncorr_errors, mult_corr_errors in nonspecial_mult:
-        # convert to from percent to fraction
-            mult_shift = (
-                1 + mult_uncorr_errors * rng.normal(size=mult_uncorr_errors.shape) / 100
-            ).prod(axis=1)
-
-            mult_shift *= (
-                1 + mult_corr_errors * rng.normal(size=(1, mult_corr_errors.shape[1])) / 100
-            ).prod(axis=1)
-
-            mult_shifts.append(mult_shift)
-            
-        #If sep_mult is true then the multiplicative shifts were not included in the covmat
-        shifts = covmat_sqrt @ rng.normal(size=covmat.shape[1])
-        mult_part = 1.
-        if sep_mult:
-            special_mult = (
-                1 + special_mult_errors * rng.normal(size=(1, 
-                special_mult_errors.shape[1])) / 100
-                ).prod(axis=1)
-            mult_part = np.concatenate(mult_shifts, axis=0)*special_mult
-        #Shifting pseudodata
-        shifted_pseudodata = (all_pseudodata + shifts)*mult_part
-        #positivity control
-        if np.all(shifted_pseudodata[full_mask] >= 0):
-            break
 
+    mult_shifts = []
+    # Prepare the per-dataset multiplicative shifts
+    for mult_uncorr_errors, mult_corr_errors in nonspecial_mult:
+    # convert to from percent to fraction
+        mult_shift = (
+            1 + mult_uncorr_errors * rng.normal(size=mult_uncorr_errors.shape) / 100
+        ).prod(axis=1)
+
+        mult_shift *= (
+            1 + mult_corr_errors * rng.normal(size=(1, mult_corr_errors.shape[1])) / 100
+        ).prod(axis=1)
+
+        mult_shifts.append(mult_shift)
+            
+    #If sep_mult is true then the multiplicative shifts were not included in the covmat
+    shifts = covmat_sqrt @ rng.normal(size=covmat.shape[1])
+    mult_part = 1.
+    if sep_mult:
+        special_mult = (
+            1 + special_mult_errors * rng.normal(size=(1, 
+        special_mult_errors.shape[1])) / 100
+        ).prod(axis=1)
+        mult_part = np.concatenate(mult_shifts, axis=0)*special_mult
+    #regularize the negative shifts
+    reg_shifts = delta(shifts, all_pseudodata, full_mask)
+    #Shifting pseudodata
+    shifted_pseudodata = (all_pseudodata + reg_shifts)*mult_part
+          
     return shifted_pseudodata
 
+def delta(shifts, data, mask):
+    new_shifts = shifts.copy()
+    neg = shifts < 0.
+    negtofix = np.logical_and(neg, mask)
+    new_shifts[negtofix] = data[negtofix] * (np.exp(shifts[negtofix] / data[negtofix]) - 1.)
+    return np.array(new_shifts)
+
+
+
 
 def indexed_make_replica(groups_index, make_replica):
     """Index the make_replica pseudodata appropriately