Skip to content

Commit 93a1ece

Browse files
mengaldomengaldo
committed
added tests
1 parent 799dbaf commit 93a1ece

File tree

2 files changed

+361
-0
lines changed

2 files changed

+361
-0
lines changed

‎tests/test_postprocessing.py‎

Lines changed: 222 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,222 @@
1+
import os
2+
import sys
3+
import shutil
4+
import subprocess
5+
import numpy as np
6+
import xarray as xr
7+
from pathlib import Path
8+
9+
# Current, parent and file paths
10+
CWD = os.getcwd()
11+
CF = os.path.realpath(__file__)
12+
CFD = os.path.dirname(CF)
13+
14+
# Import library specific modules
15+
sys.path.append(os.path.join(CFD,"../"))
16+
sys.path.append(os.path.join(CFD,"../pyspod"))
17+
from pyspod.spod_low_storage import SPOD_low_storage
18+
from pyspod.spod_low_ram import SPOD_low_ram
19+
from pyspod.spod_streaming import SPOD_streaming
20+
21+
22+
23+
def test_postprocessing_2D():
24+
25+
# data ingestion and configuration
26+
variables = ['slip_potency']
27+
file = os.path.join(CFD,'data','earthquakes_data.nc')
28+
ds = xr.open_dataset(file)
29+
t = np.array(ds['time'])
30+
x1 = np.array(ds['x'])
31+
x2 = np.array(ds['z'])
32+
X = np.array(ds[variables[0]]).T
33+
34+
# parameters
35+
params = dict()
36+
params['nt' ] = len(t)
37+
params['xdim' ] = 2
38+
params['nv' ] = 1
39+
params['dt' ] = 1
40+
params['nt' ] = t.shape[0]
41+
params['n_FFT' ] = np.ceil(32)
42+
params['n_freq' ] = params['n_FFT'] / 2 + 1
43+
params['n_overlap' ] = np.ceil(params['n_FFT'] * 50 / 100)
44+
params['savefreqs' ] = np.arange(0,params['n_freq'])
45+
params['conf_level' ] = 0.95
46+
params['n_vars' ] = 1
47+
params['n_modes_save'] = 3
48+
params['normvar' ] = False
49+
params['savedir' ] = os.path.join(CWD, 'results', Path(file).stem)
50+
params['weights' ] = np.ones([len(x1) * len(x2) * params['nv'],1])
51+
52+
# set blockwise mean
53+
params['mean'] = 'blockwise'
54+
params['savefft'] = False
55+
56+
# SPOD analysis
57+
SPOD_analysis = SPOD_low_storage(data=X, params=params, data_handler=False, variables=variables)
58+
spod = SPOD_analysis.fit()
59+
60+
# Test postprocessing and results
61+
T_approx = 12.5; tol = 1e-10
62+
freq_found, freq_idx = spod.find_nearest_freq(freq_required=1/T_approx, freq=spod.freq)
63+
modes_at_freq = spod.get_modes_at_freq(freq_idx=freq_idx)
64+
spod.plot_eigs (filename='eigs.png')
65+
spod.plot_eigs_vs_frequency(filename='eigs.png')
66+
spod.plot_eigs_vs_period (filename='eigs.png', xticks=[1, 10, 20], yticks=[1, 2, 10])
67+
spod.plot_2D_modes_at_frequency(freq_required=freq_found,
68+
freq=spod.freq,
69+
x1=x1, x2=x2,
70+
filename='modes.png')
71+
spod.plot_2D_modes_at_frequency(freq_required=freq_found,
72+
freq=spod.freq,
73+
x1=None, x2=None,
74+
equal_axes=True,
75+
filename='modes.png',
76+
plot_max=True,
77+
coastlines='regular')
78+
spod.plot_2D_modes_at_frequency(freq_required=freq_found,
79+
freq=spod.freq,
80+
x1=None, x2=None,
81+
imaginary=True,
82+
equal_axes=True,
83+
filename='modes.png',
84+
plot_max=True,
85+
coastlines='centred')
86+
spod.plot_2D_mode_slice_vs_time(freq_required=freq_found,
87+
freq=spod.freq,
88+
filename='modes.png')
89+
spod.plot_mode_tracers(freq_required=freq_found,
90+
freq=spod.freq,
91+
coords_list=[(10,10), (14,14)],
92+
filename='tracers.png')
93+
spod.plot_2D_data(time_idx=[0,10], filename='data.png')
94+
spod.plot_2D_data(time_idx=[0,10], filename='data.png', coastlines='regular')
95+
spod.plot_2D_data(time_idx=[0,10], filename='data.png', coastlines='centred')
96+
spod.plot_data_tracers(coords_list=[(10,10), (14,14)],
97+
filename='data_tracers.png')
98+
coords, idx_coords = spod.find_nearest_coords(coords=(10,10), x=[x1,x2])
99+
try:
100+
bashCmd = ["ffmpeg", " --version"]
101+
_ = subprocess.Popen(bashCmd, stdin=subprocess.PIPE)
102+
spod.generate_2D_data_video(
103+
sampling=5,
104+
time_limits=[0,t.shape[0]],
105+
filename='data_movie.mp4')
106+
spod.generate_2D_data_video(
107+
sampling=5,
108+
time_limits=[0,t.shape[0]],
109+
filename='data_movie.mp4', coastlines='regular')
110+
spod.generate_2D_data_video(
111+
sampling=5,
112+
time_limits=[0,t.shape[0]],
113+
filename='data_movie.mp4', coastlines='centred')
114+
except:
115+
print('[test_postprocessing]: ',
116+
'Skipping video making as `ffmpeg` not present.')
117+
118+
assert((np.abs(modes_at_freq[0,1,0,0]) < 8.57413617152583e-05 +tol) & \
119+
(np.abs(modes_at_freq[0,1,0,0]) > 8.57413617152583e-05 -tol))
120+
assert((np.abs(modes_at_freq[10,3,0,2]) < 0.0008816145245031309+tol) & \
121+
(np.abs(modes_at_freq[10,3,0,2]) > 0.0008816145245031309-tol))
122+
assert((np.abs(modes_at_freq[14,15,0,1]) < 0.0018284295461606808+tol) & \
123+
(np.abs(modes_at_freq[14,15,0,1]) > 0.0018284295461606808-tol))
124+
assert((np.min(np.abs(modes_at_freq)) < 8.819039169527213e-10+tol) & \
125+
(np.min(np.abs(modes_at_freq)) > 8.819039169527213e-10-tol))
126+
assert((np.max(np.abs(modes_at_freq)) < 0.28627415402845796 +tol) & \
127+
(np.max(np.abs(modes_at_freq)) > 0.28627415402845796 -tol))
128+
129+
# clean up results
130+
try:
131+
shutil.rmtree(os.path.join(CWD,'results'))
132+
except OSError as e:
133+
print("Error: %s : %s" % (os.path.join(CWD,'results'), e.strerror))
134+
try:
135+
shutil.rmtree(os.path.join(CFD,'__pycache__'))
136+
except OSError as e:
137+
print("Error: %s : %s" % (os.path.join(CFD,'__pycache__'), e.strerror))
138+
139+
140+
141+
def test_postprocessing_3D():
142+
143+
# Let's create some 2D syntetic data
144+
# and store them into a variable called p
145+
variables = ['p']
146+
x1 = np.linspace(0,10,100)
147+
x2 = np.linspace(0, 5, 50)
148+
x3 = np.linspace(0, 2, 20)
149+
xx1, xx2, xx3 = np.meshgrid(x1, x2, x3)
150+
t = np.linspace(0, 200, 1000)
151+
s_component = np.sin(xx1 * xx2 * xx3) + np.cos(xx1)**2 + np.sin(0.1*xx2) + np.sin(0.5*xx3)**2
152+
t_component = np.sin(0.1 * t)**2 + np.cos(t) * np.sin(0.5*t)
153+
p = np.empty((t_component.shape[0],)+s_component.shape)
154+
for i, t_c in enumerate(t_component):
155+
p[i] = s_component * t_c
156+
157+
# Let's define the required parameters into a dictionary
158+
params = dict()
159+
160+
# -- required parameters
161+
params['dt' ] = 1 # data time-sampling
162+
params['nt' ] = t.shape[0] # number of time snapshots (we consider all data)
163+
params['xdim' ] = 3 # number of spatial dimensions (longitude and latitude)
164+
params['nv' ] = len(variables) # number of variables
165+
params['n_FFT' ] = 100 # length of FFT blocks (100 time-snapshots)
166+
params['n_freq' ] = params['n_FFT'] / 2 + 1 # number of frequencies
167+
params['n_overlap' ] = np.ceil(params['n_FFT'] * 0 / 100) # dimension block overlap region
168+
params['mean' ] = 'blockwise' # type of mean to subtract to the data
169+
params['normalize' ] = False # normalization of weights by data variance
170+
params['savedir' ] = os.path.join(CWD, 'results', 'simple_test') # folder where to save results
171+
172+
# -- optional parameters
173+
params['weights'] = None # if set to None, no weighting (if not specified, Default is None)
174+
params['savefreqs' ] = np.arange(0,params['n_freq']) # frequencies to be saved
175+
params['n_modes_save'] = 3 # modes to be saved
176+
params['normvar' ] = False # normalize data by data variance
177+
params['conf_level' ] = 0.95 # calculate confidence level
178+
params['savefft' ] = False # save FFT blocks to reuse them in the future (saves time)
179+
180+
# Initialize libraries for the low_storage algorithm
181+
spod = SPOD_low_storage(p, params=params, data_handler=False, variables=['p'])
182+
spod.fit()
183+
184+
# Show results
185+
T_approx = 10 # approximate period = 10 days (in days)
186+
freq = spod.freq
187+
freq_found, freq_idx = spod.find_nearest_freq(freq_required=1/T_approx, freq=freq)
188+
modes_at_freq = spod.get_modes_at_freq(freq_idx=freq_idx)
189+
spod.plot_eigs (filename='eigs.png')
190+
spod.plot_eigs_vs_frequency(filename='eigs.png')
191+
spod.plot_eigs_vs_period (filename='eigs.png')
192+
spod.plot_3D_modes_slice_at_frequency(
193+
freq_required=freq_found, freq=spod.freq,
194+
x1=x1, x2=x2, x3=x3, imaginary=True, filename='modes.png', plot_max=True)
195+
spod.plot_3D_modes_slice_at_frequency(
196+
freq_required=freq_found, freq=spod.freq,
197+
x1=x1, x2=x2, x3=x3, imaginary=False,
198+
filename='modes.png', title='sim 1')
199+
spod.plot_3D_modes_slice_at_frequency(
200+
freq_required=freq_found, freq=spod.freq,
201+
x1=None, x2=None, x3=None, imaginary=False,
202+
filename='modes.png', fftshift=True,
203+
plot_max=True, equal_axes=True)
204+
spod.plot_3D_modes_slice_at_frequency(
205+
freq_required=freq_found, freq=spod.freq,
206+
x1=None, x2=None, x3=None, imaginary=False,
207+
filename='modes.png', fftshift=True,
208+
plot_max=True, slice_dim=1, equal_axes=True)
209+
spod.plot_3D_modes_slice_at_frequency(
210+
freq_required=freq_found, freq=spod.freq,
211+
x1=None, x2=None, x3=None, imaginary=True,
212+
filename='modes.png', fftshift=True,
213+
plot_max=True, slice_dim=2, equal_axes=True)
214+
spod.plot_data_tracers(coords_list=[(4,2,1)], time_limits=[0,t.shape[0]], filename='tmp.png')
215+
216+
217+
218+
219+
220+
if __name__ == "__main__":
221+
test_postprocessing_2D()
222+
test_postprocessing_3D()

‎tests/test_weights.py‎

Lines changed: 139 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,139 @@
1+
import os
2+
import sys
3+
import shutil
4+
import subprocess
5+
import numpy as np
6+
7+
# Current, parent and file paths
8+
CWD = os.getcwd()
9+
CF = os.path.realpath(__file__)
10+
CFD = os.path.dirname(CF)
11+
12+
# Import library specific modules
13+
sys.path.append(os.path.join(CFD,"../"))
14+
sys.path.append(os.path.join(CFD,"../pyspod"))
15+
from pyspod.spod_low_storage import SPOD_low_storage
16+
from pyspod.spod_low_ram import SPOD_low_ram
17+
from pyspod.spod_streaming import SPOD_streaming
18+
import pyspod.weights as weights
19+
20+
21+
22+
def test_weights_2D():
23+
24+
# Let's create some 2D syntetic data
25+
# and store them into a variable called p
26+
variables = ['p']
27+
x1 = np.linspace(0,10,100)
28+
x2 = np.linspace(0, 5, 50)
29+
xx1, xx2 = np.meshgrid(x1, x2)
30+
t = np.linspace(0, 200, 1000)
31+
s_component = np.sin(xx1 * xx2) + np.cos(xx1)**2 + np.sin(0.1*xx2)
32+
# s_component = s_component.T
33+
t_component = np.sin(0.1 * t)**2 + np.cos(t) * np.sin(0.5*t)
34+
p = np.empty((t_component.shape[0],)+s_component.shape)
35+
for i, t_c in enumerate(t_component):
36+
p[i] = s_component * t_c
37+
38+
# Let's define the required parameters into a dictionary
39+
params = dict()
40+
41+
# -- required parameters
42+
params['dt' ] = 1 # data time-sampling
43+
params['nt' ] = t.shape[0] # number of time snapshots (we consider all data)
44+
params['xdim' ] = 2 # number of spatial dimensions (longitude and latitude)
45+
params['nv' ] = len(variables) # number of variables
46+
params['n_FFT' ] = 100 # length of FFT blocks (100 time-snapshots)
47+
params['n_freq' ] = params['n_FFT'] / 2 + 1 # number of frequencies
48+
params['n_overlap' ] = np.ceil(params['n_FFT'] * 0 / 100) # dimension block overlap region
49+
params['mean' ] = 'blockwise' # type of mean to subtract to the data
50+
params['normalize' ] = True # normalization of weights by data variance
51+
params['savedir' ] = os.path.join(CWD, 'results', 'simple_test') # folder where to save results
52+
53+
# -- optional parameters
54+
params['savefreqs' ] = np.arange(0,params['n_freq']) # frequencies to be saved
55+
params['n_modes_save'] = 3 # modes to be saved
56+
params['normvar' ] = False # normalize data by data variance
57+
params['conf_level' ] = 0.95 # calculate confidence level
58+
params['savefft' ] = False # save FFT blocks to reuse them in the future (saves time)
59+
params['weights'] = weights.geo_weights_trapz_2D(lat=x2, lon=x1, R=1, n_vars=params['nv'])
60+
61+
print(p.shape)
62+
63+
# normalize data if required
64+
if params['normalize']:
65+
params['weights'] = \
66+
weights.apply_normalization(
67+
data=p,
68+
weights=params['weights'],
69+
n_variables=params['nv'],
70+
method='variance')
71+
72+
# Initialize libraries for the low_storage algorithm
73+
spod_ls = SPOD_low_storage(p, params=params, data_handler=False, variables=['p'])
74+
spod_ls.fit()
75+
76+
# Show results
77+
T_approx = 10 # approximate period = 10 days (in days)
78+
freq = spod_ls.freq
79+
freq_found, freq_idx = spod_ls.find_nearest_freq(freq_required=1/T_approx, freq=freq)
80+
modes_at_freq = spod_ls.get_modes_at_freq(freq_idx=freq_idx)
81+
tol = 1e-10
82+
assert((np.abs(modes_at_freq[5,10,0,0]) < 0.15812414887564405 +tol) & \
83+
(np.abs(modes_at_freq[5,10,0,0]) > 0.15812414887564405 -tol))
84+
assert((np.abs(modes_at_freq[0,0,0,0]) < 0.1912878813107214 +tol) & \
85+
(np.abs(modes_at_freq[0,0,0,0]) > 0.1912878813107214 -tol))
86+
assert((np.max(np.abs(modes_at_freq)) < 0.46987413376959064 +tol) & \
87+
(np.max(np.abs(modes_at_freq)) > 0.46987413376959064 -tol))
88+
89+
90+
91+
def test_weights_3D():
92+
93+
# Let's create some 2D syntetic data
94+
# and store them into a variable called p
95+
variables = ['p']
96+
x1 = np.linspace(0,10, 50)
97+
x2 = np.linspace(0, 5, 20)
98+
x3 = np.linspace(0, 2, 10)
99+
xx1, xx2, xx3 = np.meshgrid(x1, x2, x3)
100+
t = np.linspace(0, 200, 1000)
101+
s_component = np.sin(xx1 * xx2 * xx3) + np.cos(xx1)**2 + np.sin(0.1*xx2) + np.sin(0.5*xx3)**2
102+
t_component = np.sin(0.1 * t)**2 + np.cos(t) * np.sin(0.5*t)
103+
p = np.empty((t_component.shape[0],)+s_component.shape)
104+
for i, t_c in enumerate(t_component):
105+
p[i] = s_component * t_c
106+
107+
# Let's define the required parameters into a dictionary
108+
params = dict()
109+
110+
# -- required parameters
111+
params['dt' ] = 1 # data time-sampling
112+
params['nt' ] = t.shape[0] # number of time snapshots (we consider all data)
113+
params['xdim' ] = 3 # number of spatial dimensions (longitude and latitude)
114+
params['nv' ] = len(variables) # number of variables
115+
params['n_FFT' ] = 100 # length of FFT blocks (100 time-snapshots)
116+
params['n_freq' ] = params['n_FFT'] / 2 + 1 # number of frequencies
117+
params['n_overlap' ] = np.ceil(params['n_FFT'] * 0 / 100) # dimension block overlap region
118+
params['mean' ] = 'blockwise' # type of mean to subtract to the data
119+
params['normalize' ] = False # normalization of weights by data variance
120+
params['savedir' ] = os.path.join(CWD, 'results', 'simple_test') # folder where to save results
121+
122+
# -- optional parameters
123+
params['weights'] = None # if set to None, no weighting (if not specified, Default is None)
124+
params['savefreqs' ] = np.arange(0,params['n_freq']) # frequencies to be saved
125+
params['n_modes_save'] = 3 # modes to be saved
126+
params['normvar' ] = False # normalize data by data variance
127+
params['conf_level' ] = 0.95 # calculate confidence level
128+
params['savefft' ] = False # save FFT blocks to reuse them in the future (saves time)
129+
params['weights'] = weights.geo_weights_trapz_3D(lat=x2, lon=x1, R=1, z=x3, n_vars=params['nv'])
130+
131+
# Initialize libraries for the low_storage algorithm
132+
spod = SPOD_low_storage(p, params=params, data_handler=False, variables=['p'])
133+
spod.fit()
134+
135+
136+
137+
if __name__ == "__main__":
138+
test_weights_2D()
139+
test_weights_3D()

0 commit comments

Comments
 (0)