Merge 2b3a38f into 8ec128c

bin/picca_cf1d.py

@@ -61,15 +61,15 @@ def main(cmdargs):
 
     parser.add_argument('--lambda-min',
                         type=float,
-                        default=3600.,
+                        default=None,
                         required=False,
-                        help='Lower limit on observed wavelength [Angstrom]')
+                        help='Lower limit on observed wavelength [Angstrom] (default: minimum available in data)')
 
     parser.add_argument('--lambda-max',
                         type=float,
-                        default=5500.,
+                        default=None,
                         required=False,
-                        help='Upper limit on observed wavelength [Angstrom]')
+                        help='Upper limit on observed wavelength [Angstrom] (default: maximum available in data)')
 
     parser.add_argument('--z-min-sources',
                         type=float,
@@ -157,11 +157,12 @@ def main(cmdargs):
 
     # setup variables in cf
     cf.nside = args.nside
-    cf.log_lambda_min = np.log10(args.lambda_min)
-    cf.log_lambda_max = np.log10(args.lambda_max)
+    if args.lambda_min is not None:
+        cf.log_lambda_min = np.log10(args.lambda_min)
+    if args.lambda_max is not None:
+        cf.log_lambda_max = np.log10(args.lambda_max)
     cf.delta_log_lambda = args.dll
-    cf.num_pixels = int((cf.log_lambda_max - cf.log_lambda_min) /
-                        cf.delta_log_lambda + 1)
+
     cf.x_correlation = False
 
     cf.lambda_abs = constants.ABSORBER_IGM[args.lambda_abs]
@@ -221,25 +222,50 @@ def main(cmdargs):
         cf.num_data2 = num_data2
         del z_min2, z_max2
 
+
+    #here we are trying to guess reasonable lambda bounds from the data
+    if args.lambda_min is None:
+        cf.log_lambda_min=np.min([d.log_lambda[0] for hp in cf.data.values() for d in hp])
+        if args.in_dir2:
+            log_lambda_min2=np.min([d.log_lambda[0] for hp in cf.data2.values() for d in hp])
+            cf.log_lambda_min=np.min([cf.log_lambda_min,log_lambda_min2])
+        print(f"lambda_min={10**cf.log_lambda_min:.3f}")
+    if args.lambda_max is None:
+        cf.log_lambda_max=np.max([d.log_lambda[-1] for hp in cf.data.values() for d in hp])
+        if args.in_dir2:
+            log_lambda_max2=np.max([d.log_lambda[-1] for hp in cf.data2.values() for d in hp])
+            cf.log_lambda_max=np.max([cf.log_lambda_max,log_lambda_max2])
+        print(f"lambda_max={10**cf.log_lambda_max:.3f}")
+
+
+    cf.num_pixels = int((cf.log_lambda_max - cf.log_lambda_min) /
+                        cf.delta_log_lambda + 1 + 1e-10)  #+1e-7 to fix cases where the division is an integer, but numerically slightly below that int,  
+
+
+
     # Convert lists to arrays
     cf.data = {key: np.array(value) for key, value in cf.data.items()}
     if cf.x_correlation:
         cf.data2 = {key: np.array(value) for key, value in cf.data2.items()}
 
+
     # Compute the correlation function, use pool to parallelize
     cf.counter = Value('i', 0)
     cf.lock = Lock()
     context = multiprocessing.get_context('fork')
-    pool = context.Pool(processes=args.nproc)
 
     if cf.x_correlation:
         healpixs = sorted([
             key for key in list(cf.data.keys()) if key in list(cf.data2.keys())
         ])
     else:
         healpixs = sorted(list(cf.data.keys()))
-    correlation_function_data = pool.map(corr_func, healpixs)
-    pool.close()
+
+    if args.nproc>1:
+        with context.Pool(processes=args.nproc) as pool:
+            correlation_function_data = pool.map(corr_func, healpixs)
+    else:
+        correlation_function_data = [corr_func(h) for h in healpixs]
     userprint('\n')
 
     # group data from parallelisation

py/picca/cf.py

@@ -874,11 +874,13 @@ def compute_xi_1d(healpix):
                 userprint(("computing xi: {}%").format(xicounter))
             counter.value += 1
 
-        bins = ((delta.log_lambda - log_lambda_min) / delta_log_lambda +
+        select = delta.log_lambda<=log_lambda_max
+        select &= delta.log_lambda>=log_lambda_min
+        bins = ((delta.log_lambda[select] - log_lambda_min) / delta_log_lambda +
                 0.5).astype(int)
         bins = bins + num_pixels * bins[:, None]
-        delta_times_weight = delta.weights * delta.delta
-        weights = delta.weights
+        delta_times_weight = delta.weights[select] * delta.delta[select]
+        weights = delta.weights[select]
         xi1d[bins] += delta_times_weight * delta_times_weight[:, None]
         weights1d[bins] += weights * weights[:, None]
         num_pairs1d[bins] += (weights * weights[:, None] > 0.).astype(int)
@@ -912,19 +914,23 @@ def compute_xi_1d_cross(healpix):
                 userprint(("computing xi: {}%").format(xicounter))
             counter.value += 1
 
-        bins1 = ((delta1.log_lambda - log_lambda_min) / delta_log_lambda +
+        select1 = delta1.log_lambda<=log_lambda_max
+        select1 &= delta1.log_lambda>=log_lambda_min
+        bins1 = ((delta1.log_lambda[select1] - log_lambda_min) / delta_log_lambda +
                  0.5).astype(int)
-        delta_times_weight1 = delta1.weights * delta1.delta
-        weights1 = delta1.weights
+        delta_times_weight1 = delta1.weights[select1] * delta1.delta[select1]
+        weights1 = delta1.weights[select1]
 
         thingids = [delta2.thingid for delta2 in data2[healpix]]
         neighbours = data2[healpix][np.in1d(thingids, [delta1.thingid])]
         for delta2 in neighbours:
-            bins2 = ((delta2.log_lambda - log_lambda_min) / delta_log_lambda +
+            select2 = delta2.log_lambda<=log_lambda_max
+            select2 &= delta2.log_lambda>=log_lambda_min
+            bins2 = ((delta2.log_lambda[select2] - log_lambda_min) / delta_log_lambda +
                      0.5).astype(int)
             bins = bins1 + num_pixels * bins2[:, None]
-            delta_times_weight2 = delta2.weights * delta2.delta
-            weights2 = delta2.weights
+            delta_times_weight2 = delta2.weights[select2] * delta2.delta[select2]
+            weights2 = delta2.weights[select2]
             xi1d[bins] += delta_times_weight1 * delta_times_weight2[:, None]
             weights1d[bins] += weights1 * weights2[:, None]
             num_pairs1d[bins] += (weights1 * weights2[:, None] > 0.).astype(int)