Skip to content

Commit c83e619

Browse files
mengaldomengaldo
committed
solved duplicated code
1 parent 93a1ece commit c83e619

File tree

3 files changed

+116
-151
lines changed

3 files changed

+116
-151
lines changed

‎pyspod/spod_base.py‎

Lines changed: 94 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,8 @@
1111
import warnings
1212
import numpy as np
1313
import scipy.special as sc
14+
from numpy import linalg as la
15+
from scipy.fft import fft
1416

1517
# Import custom Python packages
1618
import pyspod.postprocessing as post
@@ -276,14 +278,14 @@ def n_modes(self):
276278
return self._n_modes
277279

278280
@property
279-
def n_modes_saved(self):
281+
def n_modes_save(self):
280282
'''
281283
Get the number of modes.
282284
283285
:return: the number of modes computed by the SPOD algorithm.
284286
:rtype: int
285287
'''
286-
return self._n_modes_saved
288+
return self._n_modes_save
287289

288290
@property
289291
def modes(self):
@@ -412,6 +414,11 @@ def parse_parameters(self, isrealx, params):
412414
freq[(n_DFT+1)/2+1:] = freq[(n_DFT+1)/2+1:] - (1 / dt)
413415
n_freq = len(freq)
414416

417+
n_modes_save = n_blocks
418+
if 'n_modes_save' in self._params: n_modes_save = self._params['n_modes_save']
419+
if n_modes_save > n_blocks: n_modes_save = n_blocks
420+
self._n_modes_save = n_modes_save
421+
415422
# display parameter summary
416423
print('')
417424
print('SPOD parameters')
@@ -433,6 +440,90 @@ def parse_parameters(self, isrealx, params):
433440
# ---------------------------------------------------------------------------
434441

435442

443+
# Common methods
444+
# ---------------------------------------------------------------------------
445+
446+
def compute_blocks(self, iBlk):
447+
448+
# get time index for present block
449+
offset = min(iBlk * (self._n_DFT - self._n_overlap) + self._n_DFT, self._nt) - self._n_DFT
450+
451+
# Get data
452+
Q_blk = self._data_handler(
453+
self._data, t_0=offset, t_end=self._n_DFT+offset, variables=self._variables)
454+
Q_blk = Q_blk.reshape(self._n_DFT, self._nx * self._nv)
455+
456+
# Subtract longtime or provided mean
457+
Q_blk = Q_blk[:] - self._x_mean
458+
459+
# if block mean is to be subtracted, do it now that all data is collected
460+
if self._mean_type.lower() == 'blockwise':
461+
Q_blk = Q_blk - np.mean(Q_blk, axis=0)
462+
463+
# normalize by pointwise variance
464+
if self._normvar:
465+
Q_var = np.sum((Q_blk - np.mean(Q_blk,axis=0))**2, axis=0) / (self._n_DFT-1)
466+
# address division-by-0 problem with NaNs
467+
Q_var[Q_var < 4 * np.finfo(float).eps] = 1;
468+
Q_blk = Q_blk / Q_var
469+
470+
# window and Fourier transform block
471+
self._window = self._window.reshape(self._window.shape[0],1)
472+
Q_blk = Q_blk * self._window
473+
Q_blk_hat = (self._winWeight / self._n_DFT) * fft(Q_blk, axis=0);
474+
Q_blk_hat = Q_blk_hat[0:self._n_freq,:];
475+
476+
# correct Fourier coefficients for one-sided spectrum
477+
if self._isrealx:
478+
Q_blk_hat[1:-1,:] = 2 * Q_blk_hat[1:-1,:]
479+
480+
return Q_blk_hat, offset
481+
482+
483+
484+
def compute_standard_spod(self, Q_hat_f, iFreq):
485+
486+
# compute inner product in frequency space, for given frequency
487+
M = np.matmul(Q_hat_f.conj().T, (Q_hat_f * self._weights)) / self._n_blocks
488+
489+
# extract eigenvalues and eigenvectors
490+
L,V = la.eig(M)
491+
L = np.real_if_close(L, tol=1000000)
492+
493+
# reorder eigenvalues and eigenvectors
494+
idx = np.argsort(L)[::-1]
495+
L = L[idx]
496+
V = V[:,idx]
497+
498+
# compute spatial modes for given frequency
499+
Psi = np.matmul(Q_hat_f, np.matmul(V, np.diag(1. / np.sqrt(L) / np.sqrt(self._n_blocks))))
500+
501+
# save modes in storage too in case post-processing crashes
502+
Psi = Psi[:,0:self._n_modes_save]
503+
Psi = Psi.reshape(self._xshape+(self._nv,)+(self._n_modes_save,))
504+
file_psi = os.path.join(self._save_dir_blocks,
505+
'modes1to{:04d}_freq{:04d}.npy'.format(self._n_modes_save,iFreq))
506+
np.save(file_psi, Psi)
507+
self._modes[iFreq] = file_psi
508+
self._eigs[iFreq,:] = abs(L)
509+
510+
# get and save confidence interval if required
511+
if self._conf_interval:
512+
self._eigs_c[iFreq,:,0] = self._eigs[iFreq,:] * 2 * self._n_blocks / self._xi2_lower
513+
self._eigs_c[iFreq,:,1] = self._eigs[iFreq,:] * 2 * self._n_blocks / self._xi2_upper
514+
515+
516+
517+
def store_and_save(self):
518+
self._eigs_c_u = self._eigs_c[:,:,0]
519+
self._eigs_c_l = self._eigs_c[:,:,1]
520+
file = os.path.join(self._save_dir_blocks,'spod_energy')
521+
np.savez(file, eigs=self._eigs, eigs_c_u=self._eigs_c_u, eigs_c_l=self._eigs_c_l, f=self._freq)
522+
self._n_modes = self._eigs.shape[-1]
523+
524+
# ---------------------------------------------------------------------------
525+
526+
436527

437528
# getters with arguments
438529
# ---------------------------------------------------------------------------
@@ -461,7 +552,7 @@ def get_modes_at_freq(self, freq_idx):
461552
if self._modes is None:
462553
raise ValueError('Modes not found. Consider running fit()')
463554
elif isinstance(self._modes, dict):
464-
gb_memory_modes = freq_idx * self.nx * self._n_modes_saved * \
555+
gb_memory_modes = freq_idx * self.nx * self._n_modes_save * \
465556
sys.getsizeof(complex()) * BYTE_TO_GB
466557
gb_vram_avail = psutil.virtual_memory()[1] * BYTE_TO_GB
467558
gb_sram_avail = psutil.swap_memory()[2] * BYTE_TO_GB

‎pyspod/spod_low_ram.py‎

Lines changed: 13 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,6 @@
66
import time
77
import numpy as np
88
from tqdm import tqdm
9-
from numpy import linalg as la
10-
from scipy.fft import fft
119
import shutil
1210

1311

@@ -47,57 +45,23 @@ def fit(self):
4745
# loop over number of blocks and generate Fourier realizations,
4846
# if blocks are not saved in storage
4947
if not blocks_present:
50-
5148
Q_blk = np.zeros([self._n_DFT,self._nx], dtype='complex_')
52-
5349
for iBlk in range(0,self._n_blocks):
5450

55-
# get time index for present block
56-
offset = min(iBlk * (self._n_DFT - self._n_overlap) + self._n_DFT, self._nt) - self._n_DFT
51+
# compute block
52+
Q_blk_hat, offset = self.compute_blocks(iBlk)
5753

54+
# print info file
5855
print('block '+str(iBlk+1)+'/'+str(self._n_blocks)+\
5956
' ('+str(offset)+':'+str(self._n_DFT+offset)+'); ',
6057
' Saving to directory: ', self._save_dir_blocks)
6158

62-
# # build present block
63-
# for ti in timeIdx:
64-
# x = self._X[ti]
65-
# x = x.reshape(x.size)
66-
# Q_blk[ti-offset,:] = np.subtract(x[:], self._x_mean)
67-
68-
Q_blk = self._data_handler(
69-
self._data, t_0=offset, t_end=self._n_DFT+offset, variables=self._variables)
70-
Q_blk = Q_blk.reshape(self._n_DFT, self._nx * self._nv)
71-
72-
# Subtract longtime or provided mean
73-
Q_blk = Q_blk[:] - self._x_mean
74-
75-
# if block mean is to be subtracted, do it now that all data is collected
76-
if self._mean_type.lower() == 'blockwise':
77-
Q_blk = Q_blk - np.mean(Q_blk, axis=0)
78-
79-
# normalize by pointwise variance
80-
if self._normvar:
81-
Q_var = np.sum((Q_blk - np.mean(Q_blk,axis=0))**2, axis=0) / (self._n_DFT-1)
82-
# address division-by-0 problem with NaNs
83-
Q_var[Q_var < 4 * np.finfo(float).eps] = 1;
84-
Q_blk = Q_blk / Q_var
85-
86-
# window and Fourier transform block
87-
self._window = self._window.reshape(self._window.shape[0],1)
88-
Q_blk = Q_blk * self._window
89-
Q_blk_hat = (self._winWeight / self._n_DFT) * np.fft.fft(Q_blk, axis=0);
90-
Q_blk_hat = Q_blk_hat[0:self._n_freq,:];
91-
92-
# correct Fourier coefficients for one-sided spectrum
93-
if self._isrealx:
94-
Q_blk_hat[1:-1,:] = 2 * Q_blk_hat[1:-1,:]
95-
9659
# save FFT blocks in storage memory
9760
for iFreq in range(0,self._n_freq):
9861
file = os.path.join(self._save_dir_blocks,'fft_block{:04d}_freq{:04d}.npy'.format(iBlk,iFreq))
9962
Q_blk_hat_fi = Q_blk_hat[iFreq,:]
10063
np.save(file, Q_blk_hat_fi)
64+
10165
print('------------------------------------')
10266

10367

@@ -106,19 +70,18 @@ def fit(self):
10670
print(' ')
10771
print('Calculating SPOD (low_ram)')
10872
print('------------------------------------')
109-
n_modes_save = self._n_blocks
110-
if 'n_modes_save' in self._params: n_modes_save = self._params['n_modes_save']
11173
self._eigs = np.zeros([self._n_freq,self._n_blocks], dtype='complex_')
11274
self._modes = dict()
113-
gb_memory_modes = self._n_freq * self._nx * n_modes_save * sys.getsizeof(complex()) * BYTE_TO_GB
75+
76+
gb_memory_modes = self._n_freq * self._nx * self._n_modes_save * sys.getsizeof(complex()) * BYTE_TO_GB
11477
gb_memory_avail = shutil.disk_usage(CWD)[2] * BYTE_TO_GB
11578
print('- Memory required for storing modes ~', gb_memory_modes , 'GB')
11679
print('- Available storage memory ~', gb_memory_avail , 'GB')
11780
while gb_memory_modes >= 0.99 * gb_memory_avail:
11881
print('Not enough storage memory to save all modes... halving modes to save.')
119-
n_modes_save = np.floor(n_modes_save / 2)
120-
gb_memory_modes = self._n_freq * self._nx * n_modes_save * sys.getsizeof(complex()) * BYTE_TO_GB
121-
if n_modes_save == 0:
82+
n_modes_save = np.floor(self._n_modes_save / 2)
83+
gb_memory_modes = self._n_freq * self._nx * self._n_modes_save * sys.getsizeof(complex()) * BYTE_TO_GB
84+
if self._n_modes_save == 0:
12285
raise ValueError('Memory required for storing at least one mode '
12386
'is equal or larger than available storage memory in your system ...\n'
12487
'... aborting computation...')
@@ -132,40 +95,11 @@ def fit(self):
13295
file = os.path.join(self._save_dir_blocks,'fft_block{:04d}_freq{:04d}.npy'.format(iBlk,iFreq))
13396
Q_hat_f[:,iBlk] = np.load(file)
13497

135-
# compute inner product in frequency space, for given frequency
136-
M = np.matmul(Q_hat_f.conj().T, (Q_hat_f * self._weights)) / self._n_blocks
137-
138-
# extract eigenvalues and eigenvectors
139-
L,V = la.eig(M)
140-
L = np.real_if_close(L, tol=1000000)
141-
142-
# reorder eigenvalues and eigenvectors
143-
idx = np.argsort(L)[::-1]
144-
L = L[idx]
145-
V = V[:,idx]
146-
147-
# compute spatial modes for given frequency
148-
Psi = np.matmul(Q_hat_f, np.matmul(V, np.diag(1. / np.sqrt(L) / np.sqrt(self._n_blocks))))
149-
150-
# save modes in storage too in case post-processing crashes
151-
Psi = Psi[:,0:n_modes_save]
152-
Psi = Psi.reshape(self._xshape+(self._nv,)+(n_modes_save,))
153-
file_psi = os.path.join(self._save_dir_blocks,'modes1to{:04d}_freq{:04d}.npy'.format(n_modes_save,iFreq))
154-
np.save(file_psi, Psi)
155-
self._modes[iFreq] = file_psi
156-
self._eigs[iFreq,:] = abs(L)
157-
158-
# get and save confidence interval if required
159-
if self._conf_interval:
160-
self._eigs_c[iFreq,:,0] = self._eigs[iFreq,:] * 2 * self._n_blocks / self._xi2_lower
161-
self._eigs_c[iFreq,:,1] = self._eigs[iFreq,:] * 2 * self._n_blocks / self._xi2_upper
162-
self._eigs_c_u = self._eigs_c[:,:,0]
163-
self._eigs_c_l = self._eigs_c[:,:,1]
164-
file = os.path.join(self._save_dir_blocks,'spod_energy')
165-
np.savez(file, eigs=self._eigs, eigs_c_u=self._eigs_c_u, eigs_c_l=self._eigs_c_l, f=self._freq)
166-
self._n_modes = self._eigs.shape[-1]
167-
self._n_modes_saved = n_modes_save
98+
# compute standard spod
99+
self.compute_standard_spod(Q_hat_f, iFreq)
168100

101+
# store and save results
102+
self.store_and_save()
169103

170104
# delete FFT blocks from memory if saving not required
171105
if self._savefft == False:

‎pyspod/spod_low_storage.py‎

Lines changed: 9 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,6 @@
66
import time
77
import numpy as np
88
from tqdm import tqdm
9-
from numpy import linalg as la
10-
from scipy.fft import fft
119
import psutil
1210

1311

@@ -68,37 +66,16 @@ def fit(self):
6866
Q_hat[iFreq,:,iBlk] = np.load(file)
6967
else:
7068
# loop over number of blocks and generate Fourier realizations
69+
# if blocks are not saved in storage
7170
for iBlk in range(0,self._n_blocks):
7271

73-
# get time index for present block
74-
offset = min(iBlk * (self._n_DFT - self._n_overlap) + self._n_DFT, self._nt) - self._n_DFT
72+
# compute block
73+
Q_blk_hat, offset = self.compute_blocks(iBlk)
7574

75+
# print info file
7676
print('block '+str(iBlk+1)+'/'+str(self._n_blocks)+\
7777
' ('+str(offset)+':'+str(self._n_DFT+offset)+')')
7878

79-
Q_blk = self._data_handler(
80-
self._data, t_0=offset, t_end=self._n_DFT+offset, variables=self._variables)
81-
Q_blk = Q_blk.reshape(self._n_DFT, self._nx * self._nv)
82-
Q_blk = Q_blk[:] - self._x_mean
83-
84-
# if block mean is to be subtracted, do it now that all data is collected
85-
if self._mean_type.lower() == 'blockwise':
86-
Q_blk = Q_blk - np.mean(Q_blk, axis=0)
87-
88-
# normalize by pointwise variance
89-
if self._normvar:
90-
Q_var = np.sum((Q_blk - np.mean(Q_blk,axis=0))**2, axis=0) / (self._n_DFT-1)
91-
# address division-by-0 problem with NaNs
92-
Q_var[Q_var < 4 * np.finfo(float).eps] = 1;
93-
Q_blk = Q_blk / Q_var
94-
95-
96-
# window and Fourier transform block
97-
Q_blk = Q_blk * self._window
98-
Q_blk_hat = (self._winWeight / self._n_DFT) * fft(Q_blk, axis=0)
99-
# Q_blk_hat = (winWeight / self._n_DFT) * pyfftw.interfaces.scipy_fftpack.fft(Q_blk, axis=0)
100-
Q_blk_hat = Q_blk_hat[0:self._n_freq,:];
101-
10279
# save FFT blocks in storage memory if required
10380
if self._savefft:
10481
for iFreq in range(0,self._n_freq):
@@ -107,10 +84,6 @@ def fit(self):
10784
Q_blk_hat_fi = Q_blk_hat[iFreq,:]
10885
np.save(file, Q_blk_hat_fi)
10986

110-
# correct Fourier coefficients for one-sided spectrum
111-
if self._isrealx:
112-
Q_blk_hat[1:-1,:] = 2 * Q_blk_hat[1:-1,:]
113-
11487
# store FFT blocks in RAM
11588
Q_hat[:,:,iBlk] = Q_blk_hat
11689
print('--------------------------------------')
@@ -121,10 +94,6 @@ def fit(self):
12194
print(' ')
12295
print('Calculating SPOD (low_storage)')
12396
print('--------------------------------------')
124-
n_modes_save = self._n_blocks
125-
if 'n_modes_save' in self._params: n_modes_save = self._params['n_modes_save']
126-
if n_modes_save > self._n_blocks:
127-
n_modes_save = self._n_blocks
12897
self._eigs = np.zeros([self._n_freq,self._n_blocks], dtype='complex_')
12998
self._modes = dict()
13099

@@ -134,40 +103,11 @@ def fit(self):
134103
# get FFT block from RAM memory for each given frequency
135104
Q_hat_f = np.squeeze(Q_hat[iFreq,:,:]).astype('complex_')
136105

137-
# compute inner product in frequency space, for given frequency
138-
M = np.matmul(Q_hat_f.conj().T, (Q_hat_f * self._weights)) / self._n_blocks
139-
140-
# extract eigenvalues and eigenvectors
141-
L,V = la.eig(M)
142-
L = np.real_if_close(L, tol=1000000)
143-
144-
# reorder eigenvalues and eigenvectors
145-
idx = np.argsort(L)[::-1]
146-
L = L[idx]
147-
V = V[:,idx]
148-
149-
# compute spatial modes for given frequency
150-
Psi = np.matmul(Q_hat_f, np.matmul(V, np.diag(1. / np.sqrt(L) / np.sqrt(self._n_blocks))))
151-
152-
# save modes in storage too in case post-processing crashes
153-
Psi = Psi[:,0:n_modes_save]
154-
Psi = Psi.reshape(self._xshape+(self._nv,)+(n_modes_save,))
155-
file_psi = os.path.join(self._save_dir_blocks,
156-
'modes1to{:04d}_freq{:04d}.npy'.format(n_modes_save,iFreq))
157-
np.save(file_psi, Psi)
158-
self._modes[iFreq] = file_psi
159-
self._eigs[iFreq,:] = abs(L)
160-
161-
# get and save confidence interval if required
162-
if self._conf_interval:
163-
self._eigs_c[iFreq,:,0] = self._eigs[iFreq,:] * 2 * self._n_blocks / self._xi2_lower
164-
self._eigs_c[iFreq,:,1] = self._eigs[iFreq,:] * 2 * self._n_blocks / self._xi2_upper
165-
self._eigs_c_u = self._eigs_c[:,:,0]
166-
self._eigs_c_l = self._eigs_c[:,:,1]
167-
file = os.path.join(self._save_dir_blocks,'spod_energy')
168-
np.savez(file, eigs=self._eigs, eigs_c_u=self._eigs_c_u, eigs_c_l=self._eigs_c_l, f=self._freq)
169-
self._n_modes = self._eigs.shape[-1]
170-
self._n_modes_saved = n_modes_save
106+
# compute standard spod
107+
self.compute_standard_spod(Q_hat_f, iFreq)
108+
109+
# store and save results
110+
self.store_and_save()
171111
print('--------------------------------------')
172112
print(' ')
173113

0 commit comments

Comments
 (0)