updated arc_sic with faster interpolation

algorithm_config.yaml

 description: Estimates canopy height from InSAR coherence and LiDAR data
 algo_name: ich
-version: 0.4
+version: 0.4.1
 environment: ubuntu
 repository_url: https://repo.maap-project.org/bnarayanarao/insar_forest_height.git
 docker_url: mas.maap-project.org/root/maap-workspaces/base_images/python:v4.1.0

src/ich/algo.py

@@ -24,95 +24,86 @@ from scipy.interpolate import interpn
 from sklearn.metrics import mean_squared_error as mse
 from concurrent.futures import ProcessPoolExecutor, as_completed
 import os
+from concurrent.futures import ThreadPoolExecutor
+from scipy.interpolate import PchipInterpolator
+
+# Precompute XX and YY
+XX = np.linspace(0, np.pi, num=100, endpoint=True)
+XX[0] = np.spacing(1)  # Avoid division by zero
+YY = np.sin(XX) / XX
+YY[0] = 1  # Correct first value
+YY[-1] = 0  # Correct last value
+XX = XX[::-1]  # Flip arrays
+YY = YY[::-1]
+
+# Precompute the interpolation function
+interp_func = PchipInterpolator(YY, XX)
+
+def arc_sinc_fast(x, c_param):
+    # Clip x to the range [0, 1], works for both 1D and 2D arrays
+    x_clipped = np.clip(x, 0, 1)
+    
+    # Interpolate and apply c_param scaling (works for both 1D and 2D arrays)
+    y = interp_func(x_clipped) * c_param
+    
+    # Clip final result to the range [0, inf), works for both 1D and 2D arrays
+    return np.maximum(y, 0)
 
-def arc_sinc(x, c_param):
-    x = np.clip(x, 0, 1)  # Ensure x is between 0 and 1
+
+def arc_sinc_(x, c_param):
+    # Get rid of extreme values by set all values where x > 1 equal to 1, and x < 0 equal to 0 
+    if x >1:
+        x=1
+    elif x<0:
+        x=0
+    # x[(x > 1)] = 1
+    # x[(x < 0)] = 0
+
+    # Create array of increments between 0 and pi of size pi/100
     XX = np.linspace(0, np.pi, num=100, endpoint=True)
-    XX[0] = np.spacing(1)  # Avoid division by zero issues
+
+    # Set the first value of XX to eps to avoid division by zero issues -> Paul's suggestion
+    XX[0] = np.spacing(1)
+
+    # Calculate sinc for XX and save it to YY
+    ## YY = sinc(XX / math.pi)
     YY = np.sin(XX) / XX
+
+    # Reset the first value of XX to zero and the first value of YY to the corresponding output
     XX[0] = 0
     YY[0] = 1
+    
+    # Set the last value of YY to 0 to avoid NaN issues
     YY[-1] = 0
+    # Flip XX and YY left to right
     XX = XX[::-1]
     YY = YY[::-1]
-    interp_func = interpolate.interp1d(YY, XX * c_param, kind='slinear', bounds_error=False, fill_value=0)
+    # print(XX.shape,YY.shape)
+     # Run interpolation
+    # XX and YY are your original values, x is the query values, and y is the interpolated values that correspond to x
+    interp_func = interpolate.interp1d(YY, XX * c_param, kind='slinear') 
     y = interp_func(x)
-    return np.clip(y, 0, None)  # Ensure no negative heights
-
-from concurrent.futures import ThreadPoolExecutor
-import os
-
-
-def arc_sinc_block(x_chunk, c_param_chunk, XX, YY):
-    # Initialize result array for the chunk
-    y_chunk = np.zeros_like(x_chunk)
-    
-    # Handle the interpolation for each unique c_param value in the chunk
-    for unique_c in np.unique(c_param_chunk):
-        # Interpolation function for the current unique c_param
-        interp_func = interpolate.interp1d(YY, XX * unique_c, kind='slinear', bounds_error=False, fill_value=0)
-        
-        # Apply the interpolation
-        mask = (c_param_chunk == unique_c)
-        y_chunk[mask] = interp_func(x_chunk[mask])
-    
-    # Clip the result to ensure non-negative values
-    return np.clip(y_chunk, 0, None)
-
 
+    # Set all values in y less than 0 equal to 0
+    y[(y < 0)] = 0
+    # return y
+    return y
 
 
-def arc_sinc_vectorized_parallel(x, c_param, chunk_size=100000):
-    # Ensure x and c_param are numpy arrays
-    x = np.asarray(x)
-    c_param = np.asarray(c_param)
-    
-    # Clip x to the range [0, 1]
-    x = np.clip(x, 0, 1)
-    
-    # Create array of increments between 0 and pi
+def arc_sinc(x, c_param):
+    x = np.clip(x, 0, 1)  # Ensure x is between 0 and 1
     XX = np.linspace(0, np.pi, num=100, endpoint=True)
-    
-    # Set the first value of XX to eps to avoid division by zero issues
-    XX[0] = np.spacing(1)
-    
-    # Calculate sinc for XX and save it to YY
+    XX[0] = np.spacing(1)  # Avoid division by zero issues
     YY = np.sin(XX) / XX
-    
-    # Reset the first value of XX and YY
     XX[0] = 0
     YY[0] = 1
     YY[-1] = 0
-    
-    # Reverse arrays for interpolation
     XX = XX[::-1]
     YY = YY[::-1]
-    
-    # Get the length of the input arrays
-    n = len(x)
-    
-    # Automatically get the number of available CPU cores
-    num_workers = os.cpu_count()
-    
-    # Split input arrays into chunks and process in parallel
-    y = np.zeros_like(x)
-    with ThreadPoolExecutor(max_workers=num_workers) as executor:
-        futures = []
-        
-        # Loop over chunks
-        for i in range(0, n, chunk_size):
-            x_chunk = x[i:i + chunk_size]
-            c_param_chunk = c_param[i:i + chunk_size]
-            
-            # Submit chunk processing to the thread pool
-            futures.append(executor.submit(arc_sinc_block, x_chunk, c_param_chunk, XX, YY))
-        
-        # Collect results
-        for idx, future in enumerate(futures):
-            i = idx * chunk_size
-            y[i:i + chunk_size] = future.result()
-    
-    return y
+    interp_func = interpolate.interp1d(YY, XX * c_param, kind='slinear', bounds_error=False, fill_value=0)
+    y = interp_func(x)
+    return np.clip(y, 0, None)  # Ensure no negative heights
+
 
 def arc_sinc_vectorized(x, c_param):
     # Ensure x and c_param are numpy arrays
@@ -170,6 +161,9 @@ def cal_(cor, gedi, htl, htg):
             for C_param in C_range:
                 gama = np.clip(temp_cor / S_param, 0, 1)
                 height = arc_sinc(gama, C_param)
+                # height = arc_sinc_fast(gama, C_param)
+
+                
                 height[(height >= htg) | (height <= htl)] = np.nan
                 x, y = temp_gedi.flatten(), height.flatten()
                 valid = ~np.isnan(x) & ~np.isnan(y)
@@ -204,6 +198,7 @@ def cal_(cor, gedi, htl, htg):
                 for C_param in np.arange(cCoarse - 1, cCoarse + 1, 0.2):
                     gama = np.clip(temp_cor / S_param, 0, 1)
                     height = arc_sinc(gama, C_param)
+                    # height = arc_sinc_fast(gama, C_param)
                     height[(height >= htg) | (height <= htl)] = np.nan
                     x, y = temp_gedi.flatten(), height.flatten()
                     valid = ~np.isnan(x) & ~np.isnan(y)
@@ -267,6 +262,7 @@ def process_block(i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_):
 
             gama = cohBlock / parm[1]
             ht = arc_sinc(gama, parm[2]) * mask
+            # ht = arc_sinc_fast(gama, parm[2]) * mask
 
             return start_i, end_i, start_j, end_j, s_parm, c_parm, rmse_parm, ht, count
         else:

src/ich/args_in.py

@@ -26,7 +26,7 @@ from scipy.interpolate import interpn
 from empatches import EMPatches
 from sklearn.model_selection import train_test_split
 from plotting import density_scatter, plot_data, rmse
-from algo import arc_sinc, arc_sinc_vectorized, cal_, dynamicWindow,arc_sinc_vectorized_parallel
+from algo import arc_sinc_fast,arc_sinc, arc_sinc_, cal_, dynamicWindow
 from rst_io import read_bin, write_tif, write_bin
 from utils import blockshaped, unblockshaped, unpad, spatial_intp_lin, split_sample
 
@@ -191,7 +191,14 @@ def rvog_inverse(args):
     height = []
     print('...')
     # height = arc_sinc_vectorized(gama, c_temp)
-    height = arc_sinc_vectorized_parallel(gama, c_temp)
+    # height = arc_sinc_vectorized_parallel(gama, c_temp)
+    height = arc_sinc_fast(gama, c_temp)
+
+    # for ii in tqdm(range(np.size(c_temp))):
+    # # for ii in (range(np.size(c_temp))):
+    #     height.append(arc_sinc_(gama[ii], c_temp[ii]))
+
+
     height = np.array(height).reshape(np.shape(temp_cor))
     height[height < 0] = 0
     print("%0.2f sec"%(time.time()-t2))