From af5619c91d10d3ad510930a93a125bdc845d3b51 Mon Sep 17 00:00:00 2001
From: Narayana-Rao <narayanarao.bhogapurapu@gmail.com>
Date: Mon, 2 Dec 2024 18:31:33 -0800
Subject: [PATCH] memory issue fix
---
src/sarbio/algo.py | 123 ++++++++++++++++++++++++++++++++++++------
src/sarbio/args_in.py | 62 +++++++++++++--------
2 files changed, 148 insertions(+), 37 deletions(-)
diff --git a/src/sarbio/algo.py b/src/sarbio/algo.py
index 76e6a40..6da28cd 100644
--- a/src/sarbio/algo.py
+++ b/src/sarbio/algo.py
@@ -127,6 +127,88 @@ def process_cal_window(win, temp_sig, temp_inc, temp_er, temp_s, temp_lidar, li
#################################################
##################################################
+def inversion__temp(temp_sig, temp_inc, temp_er, temp_s, temp_agbd, A, B, C, alpha, beta, gama, agbdl, agbdg):
+ # Generate biomass step range
+ biom_step = np.arange(agbdl, agbdg + 0.1, 0.1).astype(np.float32)
+
+
+ dim = np.shape(temp_inc)
+ k = 2 * np.pi / 0.23 * np.ones_like(temp_inc.flatten())
+ # slices = 10
+
+ slices = max(10,int(np.size(k)//1e5))
+ # Handle padding
+ rPad = len(temp_inc.flatten()) % slices
+ if rPad != 0:
+ padding = slices - rPad
+ pad_mode = 'constant'
+ pad_val = np.nan
+ inc_ = np.pad(temp_inc.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ sigma_vh_ = np.pad(temp_sig.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ er_ = np.pad(temp_er.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ s_ = np.pad(temp_s.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ k_ = np.pad(k.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ A_ = np.pad(A.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ B_ = np.pad(B.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ C_ = np.pad(C.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ alpha_ = np.pad(alpha.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ beta_ = np.pad(beta.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ gama_ = np.pad(gama.flatten(), (0, padding), mode=pad_mode, constant_values=pad_val)
+ else:
+ inc_, sigma_vh_, er_, s_, k_, A_, B_, C_, alpha_, beta_, gama_ = map(lambda arr: arr.flatten(), [temp_inc, temp_sig, temp_er, temp_s, k,A, B, C, alpha, beta, gama])
+
+ # Split arrays into slices
+ inc_, sigma_vh_, er_, s_, k_, A_, B_, C_, alpha_, beta_, gama_ = map(lambda arr: np.array_split(arr, slices),
+ [inc_, sigma_vh_, er_, s_, k_, A_, B_, C_, alpha_, beta_, gama_])
+
+ biom_est = []
+ error_list = []
+ est_sig = []
+
+ # Loop to compute biomass estimation
+ for i in range(np.shape(inc_)[0]):
+ # Prepare variables for each slice
+ inc_in, sigvh_in, er_in, k_in, s_in, A_in, B_in, C_in, alpha_in, beta_in, gama_in = map(
+ lambda x: np.tile(x[i], (np.size(biom_step), 1)).T,
+ [inc_, sigma_vh_, er_, k_, s_, A_, B_, C_, alpha_, beta_, gama_])
+
+ # Biomass simulation
+ bio_sim = np.tile(biom_step, (np.size(inc_[i].flatten()), 1))
+ sigvh_sim = NISAR_bio_HV((bio_sim, er_in, inc_in, k_in, s_in), A_in, B_in, C_in, alpha_in, beta_in, gama_in)
+
+ # Calculate error
+ error = np.abs(sigvh_in - sigvh_sim)
+ ind1 = np.arange(len(bio_sim))
+ ind2 = np.argmin(error, axis=1)
+ print("error_shape",error.shape)
+ print("sigvh_in_shape",sigvh_in.shape)
+ print("sigvh_sim_shape",sigvh_sim.shape)
+ # Collect results
+ biom_est.append(bio_sim[ind1, ind2])
+ error_list.append(error[ind1, ind2])
+ est_sig.append(sigvh_sim[ind1, ind2])
+ print(bio_sim[ind1, ind2].shape, error[ind1, ind2].shape, sigvh_sim[ind1, ind2].shape)
+
+ # Concatenate and reshape the results
+ biom_est = np.concatenate(biom_est)
+ sigma_vh = np.concatenate(sigma_vh_)
+ error_list = np.concatenate(error_list)
+ est_sig = np.concatenate(est_sig)
+
+ # Adjust length if necessary
+ if rPad > 0:
+ biom_est = biom_est[:-(slices - rPad)]
+ sigma_vh = sigma_vh[:-(slices - rPad)]
+ error_list = error_list[:-(slices - rPad)]
+ est_sig = est_sig[:-(slices - rPad)]
+
+ # Reshape arrays to match the original dimensions
+ biom_est = biom_est.reshape(dim)
+ error_list = error_list.reshape(dim)
+ est_sig = est_sig.reshape(dim)
+
+ return biom_est, error_list, est_sig
+
def inversion(temp_sig, temp_inc, temp_er, temp_s, temp_agbd, A,B,C,alpha,beta,gama,agbdl, agbdg):
biom_step = np.arange(agbdl,agbdg+.1,.1).astype(np.float32)
@@ -134,11 +216,10 @@ def inversion(temp_sig, temp_inc, temp_er, temp_s, temp_agbd, A,B,C,alpha,beta,g
dim = np.shape(temp_inc)
k = np.repeat(2*np.pi/.23,np.size(temp_inc.flatten()))
- slices = 10
+ # slices = 10
+ slices = max(10,int(np.size(k)//25e3))
rPad = len(temp_inc.flatten()) %slices
- # print('pad:', rPad)
- k = np.repeat(2*np.pi/.23,np.size(temp_inc.flatten()))
-
+
# print(np.nanmean(temp_sig),np.nanmean(sigma_model))
if rPad != 0:
inc_ = np.pad(temp_inc.flatten(), (0, slices - rPad), mode='constant', constant_values=np.nan)
@@ -165,7 +246,7 @@ def inversion(temp_sig, temp_inc, temp_er, temp_s, temp_agbd, A,B,C,alpha,beta,g
alpha_ = alpha.flatten().copy()
beta_ = beta.flatten().copy()
gama_ = gama.flatten().copy()
-
+ del temp_inc, temp_sig, temp_er, temp_s, temp_agbd, A,B,C,alpha,beta,gama,k
inc_ = np.array_split(inc_, slices)
sigma_vh_ = np.array_split(sigma_vh_, slices)
er_ = np.array_split(er_, slices)
@@ -181,6 +262,7 @@ def inversion(temp_sig, temp_inc, temp_er, temp_s, temp_agbd, A,B,C,alpha,beta,g
biom_est=[]
error_list=[]
est_sig = []
+ # print('inc_shape',np.shape(inc_))
for i in range(np.shape(inc_)[0]):
inc_in = np.tile(inc_[i], (np.size(biom_step), 1)).T
sigvh_in = np.tile(sigma_vh_[i], (np.size(biom_step), 1)).T
@@ -197,14 +279,24 @@ def inversion(temp_sig, temp_inc, temp_er, temp_s, temp_agbd, A,B,C,alpha,beta,g
bio_sim = np.tile(biom_step, (np.size(inc_[i].flatten()), 1))
sigvh_sim = NISAR_bio_HV((bio_sim,er_in,inc_in,k_in,s_in),A_in,B_in,C_in,alpha_in,beta_in,gama_in)
- error = np.abs(sigvh_in-sigvh_sim)
- ind1 = np.arange(len(bio_sim))
- ind2 = np.argmin(error, axis=1)
- result = bio_sim[ind1, ind2]
- biom_est.append(result)
- error_list.append(error[ind1, ind2])
- est_sig.append(sigvh_sim[ind1, ind2])
-
+ try:
+ error = np.abs(sigvh_in-sigvh_sim)
+ # print("error_shape",error.shape)
+ # print("sigvh_in_shape",sigvh_in.shape)
+ # print("sigvh_sim_shape",sigvh_in.shape)
+ ind1 = np.arange(len(bio_sim))
+ ind2 = np.argmin(error, axis=1)
+
+ biom_est.append(bio_sim[ind1, ind2])
+ error_list.append(error[ind1, ind2])
+ est_sig.append(sigvh_sim[ind1, ind2])
+ # print("bio_sim",bio_sim[ind1, ind2].shape, "error",error[ind1, ind2].shape, "sigvh_sim",sigvh_sim[ind1, ind2].shape)
+ # print('inc_shape',np.shape(inc_))
+ except:
+ biom_est.append(np.repeat(np.nan,np.shape(inc_)[1]))
+ error_list.append(np.repeat(np.nan,np.shape(inc_)[1]))
+ est_sig.append(np.repeat(np.nan,np.shape(inc_)[1]))
+
biom_est = np.concatenate(biom_est)
if rPad>0:
biom_est = biom_est[:-(slices-rPad)]
@@ -284,12 +376,13 @@ def calibrate(temp_sig, temp_inc, temp_er, temp_s, temp_agbd, agbdl, agbdg):
bounds=((0,0,0,0,0,0),(1,1,1,1,1,1)),\
method='dogbox', maxfev=100000,\
)
- except:
+ except Exception as e:
+ # print(f"Exception in cal_: {e}")
nlfit = np.array([0.013,0.00195,0.0047, 0.11,0.9,0.23])
if np.all(nlfit == 0):
nlfit[:] = np.nan
- # # print(nlfit)
+ # print(nlfit)
# sigma_model = NISAR_bio_HV((biomass,er,inc,k,s),nlfit[0],nlfit[1],nlfit[2],nlfit[3],nlfit[4],nlfit[5])
# try:
diff --git a/src/sarbio/args_in.py b/src/sarbio/args_in.py
index 6a7c0e6..d51a867 100644
--- a/src/sarbio/args_in.py
+++ b/src/sarbio/args_in.py
@@ -64,7 +64,8 @@ def agb_inverse(args):
))
- lidar_agbd = read_bin(args.lidarFile)
+ lidar_agbd = read_bin(args.lidarFile)
+ lidar_agbd[lidar_agbd<=0] = np.nan
backscatterFile = args.backscatterFile
sigma_hv = 10**(read_bin(backscatterFile)/10)
inc = read_bin(args.incFile)/100
@@ -105,13 +106,13 @@ def agb_inverse(args):
pp99 = np.nanpercentile(lidar_agbd,99.9)
args.agbdg = int( 5 * np.ceil( pp99 / 5. ))
-
lidar_agbd[lidar_agbd<=args.agbdl] = np.nan
lidar_agbd[lidar_agbd>args.agbdg] = np.nan
lidar_Nsamples = np.sum(~np.isnan(lidar_agbd))
# sigma_hv[sigma_hv>0]=np.nan
rows,cols = np.shape(sigma_hv)[0],np.shape(sigma_hv)[1]
+ print(f"lidar_Nsamples: {lidar_Nsamples} agbdl: {args.agbdl} agbdg: {args.agbdg}")
######################################################
if args.validation>0 and args.validation<1:
@@ -219,7 +220,7 @@ def agb_inverse(args):
# Unpack results
A, B,C,alpha,beta,gama, count = zip(*results)
-
+ del results
temp_mask = np.zeros(temp_lidar.shape)
for win in range(np.shape(temp_lidar)[0]):
mask = temp_lidar[win,:,:].copy()
@@ -284,8 +285,23 @@ def agb_inverse(args):
gamai = spatial_intp_lin(gama)
print("%0.2f sec "%(time.time()-t2))
-
-
+ t3 = time.time()
+ if rPad!=0 or cPad!=0:
+ A = unpad(A, pad_width)
+ B = unpad(B, pad_width)
+ C = unpad(C, pad_width)
+ alpha = unpad(alpha, pad_width)
+ beta = unpad(beta, pad_width)
+ gama = unpad(gama, pad_width)
+ write_tif(outAFile,A,os.path.basename(args.backscatterFile))
+ write_tif(outBFile,B,os.path.basename(args.backscatterFile))
+ write_tif(outCFile,C,os.path.basename(args.backscatterFile))
+
+ write_tif(outalphaFile,alpha,os.path.basename(args.backscatterFile))
+ write_tif(outbetaFile,beta,os.path.basename(args.backscatterFile))
+ write_tif(outgamaFile,gama,os.path.basename(args.backscatterFile))
+ del A,B,C,alpha,beta,gama
+
Ai = blockshaped(Ai, args.window_size, args.window_size)
Bi = blockshaped(Bi, args.window_size, args.window_size)
Ci = blockshaped(Ci, args.window_size, args.window_size)
@@ -301,7 +317,7 @@ def agb_inverse(args):
batch_size = max(num_windows//16,100)
results = [None] * num_windows # Preallocate a list to hold results
- with concurrent.futures.ProcessPoolExecutor() as executor:
+ with concurrent.futures.ProcessPoolExecutor(max_workers=os.cpu_count()-1) as executor:
futures = {}
for start in range(0, num_windows, batch_size):
end = min(start + batch_size, num_windows)
@@ -319,7 +335,9 @@ def agb_inverse(args):
biom_est, error_list, est_sig = zip(*results)
-
+ del results,temp_mask,temp_sig,temp_inc, temp_er, temp_s,temp_lidar
+
+
Ai = unblockshaped(Ai, rows,cols)
Bi = unblockshaped(Bi, rows,cols)
Ci = unblockshaped(Ci, rows,cols)
@@ -550,18 +568,18 @@ def agb_inverse(args):
if rPad!=0 or cPad!=0:
lidar_agbd_cal = unpad(lidar_agbd_cal, pad_width)
- temp_sig = unpad(temp_sig, pad_width)
- A = unpad(A, pad_width)
- B = unpad(B, pad_width)
- C = unpad(C, pad_width)
- alpha = unpad(alpha, pad_width)
- beta = unpad(beta, pad_width)
- gama = unpad(gama, pad_width)
+ # temp_sig = unpad(temp_sig, pad_width)
+ # A = unpad(A, pad_width)
+ # B = unpad(B, pad_width)
+ # C = unpad(C, pad_width)
+ # alpha = unpad(alpha, pad_width)
+ # beta = unpad(beta, pad_width)
+ # gama = unpad(gama, pad_width)
# temp_sig = unpad(temp_sig, pad_width)
# temp_er = unpad(temp_er, pad_width)
# temp_inc = unpad(temp_inc, pad_width)
# temp_s = unpad(temp_s, pad_width)
- temp_lidar = unpad(temp_lidar, pad_width)
+ # temp_lidar = unpad(temp_lidar, pad_width)
Ai = unpad(Ai, pad_width)
Bi = unpad(Bi, pad_width)
Ci = unpad(Ci, pad_width)
@@ -575,7 +593,7 @@ def agb_inverse(args):
# print("%0.2f sec"%(time.time()-t1))
- del temp_mask,temp_sig,temp_inc, temp_er, temp_s
+
plt_agbd_est = biom_est.copy()
plt_agbd_est[plt_agbd_est<=1] = np.nan
plt_agbd_est[plt_agbd_est>args.agbdg-5] = np.nan
@@ -595,13 +613,13 @@ def agb_inverse(args):
print("[4/4] Writing output files... ",end=" ")
- write_tif(outAFile,A,os.path.basename(args.backscatterFile))
- write_tif(outBFile,B,os.path.basename(args.backscatterFile))
- write_tif(outCFile,C,os.path.basename(args.backscatterFile))
+ # write_tif(outAFile,A,os.path.basename(args.backscatterFile))
+ # write_tif(outBFile,B,os.path.basename(args.backscatterFile))
+ # write_tif(outCFile,C,os.path.basename(args.backscatterFile))
- write_tif(outalphaFile,alpha,os.path.basename(args.backscatterFile))
- write_tif(outbetaFile,beta,os.path.basename(args.backscatterFile))
- write_tif(outgamaFile,gama,os.path.basename(args.backscatterFile))
+ # write_tif(outalphaFile,alpha,os.path.basename(args.backscatterFile))
+ # write_tif(outbetaFile,beta,os.path.basename(args.backscatterFile))
+ # write_tif(outgamaFile,gama,os.path.basename(args.backscatterFile))
write_tif(outAiFile,Ai,os.path.basename(args.backscatterFile))
write_tif(outBiFile,Bi,os.path.basename(args.backscatterFile))
--
GitLab