Skip to content

Commit

Permalink
Test all examples and fix various small bugs.
Browse files Browse the repository at this point in the history
  • Loading branch information
matthew-gaddes committed Nov 30, 2021
1 parent a4787e5 commit ea0e217
Show file tree
Hide file tree
Showing 58 changed files with 61 additions and 14 deletions.
61 changes: 53 additions & 8 deletions example_spatial_01.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,40 @@
from icasar.icasar_funcs import ICASAR
from icasar.aux import col_to_ma, r2_to_r3

def generate_random_ifg_dates(n_dates):
"""A function to generate random synthetic interferogram dates
"""

def daisy_chain_from_acquisitions(acquisitions):
"""Given a list of acquisiton dates, form the names of the interferograms that would create a simple daisy chain of ifgs.
Inputs:
acquisitions | list | list of acquistiion dates in form YYYYMMDD
Returns:
daisy_chain | list | names of daisy chain ifgs, in form YYYYMMDD_YYYYMMDD
History:
2020/02/16 | MEG | Written
"""
daisy_chain = []
n_acqs = len(acquisitions)
for i in range(n_acqs-1):
daisy_chain.append(f"{acquisitions[i]}_{acquisitions[i+1]}")
return daisy_chain

#import pdb; pdb.st_trace()

import datetime as dt
import random
day0 = "20150101"
day0_dt = dt.datetime.strptime(day0, '%Y%m%d')
acq_dates_dt = [day0_dt]
for i in range(n_dates):
t_baseline = random.choice([6,12,18,24])
acq_dates_dt.append(acq_dates_dt[-1] + dt.timedelta(days = t_baseline))
acq_dates = [dt.datetime.strftime(i, '%Y%m%d') for i in acq_dates_dt]

ifg_dates = daisy_chain_from_acquisitions(acq_dates)
return ifg_dates


#%% Things to set

Expand All @@ -27,6 +61,7 @@
"ica_param" : (1e-2, 150), # (tolerance, max iterations)
"hdbscan_param" : (35,10), # (min_cluster_size, min_samples) Discussed in more detail in Mcinnes et al. (2017). min_cluster_size sets the smallest collection of points that can be considered a cluster. min_samples sets how conservative the clustering is. With larger values, more points will be considered noise.
"out_folder" : Path('example_spatial_01_outputs'), # outputs will be saved here
"load_fastICA_results" : True,
"figures" : "png+window"} # if png, saved in a folder as .png. If window, open as interactive matplotlib figures,
# if 'png+window', both.
# default is "window" as 03_clustering_and_manifold is interactive.
Expand Down Expand Up @@ -54,14 +89,14 @@
axes[1,i].plot(range(A_dc.shape[0]), A_dc[:,i])
axes[1,i].axhline(0)
fig1.suptitle('Synthetic sources and time courses')
fig1.canvas.set_window_title("Synthetic sources and time courses")
fig1.canvas.manager.set_window_title("Synthetic sources and time courses")


fig2, axes = plt.subplots(2,5) # plot the synthetic interferograms
for i, ax in enumerate(np.ravel(axes[:])):
ax.imshow(col_to_ma(phUnw[i,:], pixel_mask))
fig2.suptitle('Mixtures (intererograms)')
fig2.canvas.set_window_title("Mixtures (intererograms)")
fig2.canvas.manager.set_window_title("Mixtures (intererograms)")

fig3, axes = plt.subplots(1,3, figsize = (11,4)) # plot a schematic of how the data are organised
axes[0].imshow(X_dc, aspect = 500)
Expand All @@ -70,19 +105,28 @@
axes[1].set_title('Mask')
axes[2].imshow(col_to_ma(X_dc[0,:], pixel_mask))
axes[2].set_title('Interferogram 1')
fig3.canvas.set_window_title("Interferograms as row vectors and a mask")
fig3.canvas.manager.set_window_title("Interferograms as row vectors and a mask")

#%% do ICA with ICSAR function

spatial_data = {'mixtures_r2' : phUnw,
# note that ICASAR now requires the ifg_dates to be supplied, so we must synthetise these too.

ifg_dates_dc = generate_random_ifg_dates(phUnw.shape[0])

spatial_data = {'ifgs_dc' : phUnw,
'mask' : pixel_mask,
'lons' : lons, # for the simplest case, these aren't needed
'lats' : lats} # for the simplest case, these aren't needed
'lats' : lats, # for the simplest case, these aren't needed
'ifg_dates_dc' : ifg_dates_dc}




S_best, time_courses, x_train_residual_ts, Iq, n_clusters, S_all_info, phUnw_mean = ICASAR(spatial_data = spatial_data, **ICASAR_settings)


#%% We can reconstruct the data using the sources and timecourses, but don't forget that ICA returns mean centered sources
phUnw_mean = phUnw_mean[:, np.newaxis]

X_dc_reconstructed = (time_courses @ S_best) + phUnw_mean # here we add the mean back
X_dc_reconstructed_source0 = (time_courses[:,0:1] @ S_best[0:1,:]) + phUnw_mean # and remake the entire time series using only IC0
Expand All @@ -101,18 +145,19 @@
axes[2].set_title('Reconstructed Ifg. \n (IC0 only)')
fig4.colorbar(im2, ax = axes[2])

fig4.canvas.set_window_title("Reconstructed Data")
fig4.canvas.manager.set_window_title("Reconstructed Data")


#%% Note that the amount of bootstrapping done by ICASAR can also be controlled, and seen in the clustering and 2d manifold plot:

ICASAR_settings["bootstrapping_param"] = (100, 100) # (number of runs with bootstrapping, number of runs without bootstrapping)
ICASAR_settings['out_folder'] = 'example_spatial_01_outputs_part2'

spatial_data = {'mixtures_r2' : phUnw,
spatial_data = {'ifgs_dc' : phUnw,
'mask' : pixel_mask,
'lons' : lons,
'lats' : lats}
'lats' : lats,
'ifg_dates_dc': ifg_dates_dc}

S_best, time_courses, x_train_residual_ts, Iq, n_clusters, S_all_info, phUnw_mean = ICASAR(spatial_data = spatial_data, **ICASAR_settings)

Binary file modified example_spatial_01_outputs/01_pca_variance_line.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
Binary file not shown.
Binary file modified example_spatial_01_outputs/ICASAR_results.pkl
Binary file not shown.
Binary file modified example_spatial_01_outputs/ICs.kmz
Binary file not shown.
Binary file modified example_spatial_01_outputs_part2/01_pca_variance_line.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
Binary file not shown.
Binary file modified example_spatial_01_outputs_part2/ICASAR_results.pkl
Binary file not shown.
Binary file modified example_spatial_01_outputs_part2/ICs.kmz
Binary file not shown.
4 changes: 2 additions & 2 deletions example_spatial_02.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@
"ica_param" : (1e-2, 150), # (tolerance, max iterations)
"hdbscan_param" : (100,10), # (min_cluster_size, min_samples) Discussed in more detail in Mcinnes et al. (2017). min_cluster_size sets the smallest collection of points that can be considered a cluster. min_samples sets how conservative the clustering is. With larger values, more points will be considered noise.
"out_folder" : Path('example_spatial_02_outputs_sICA'), # outputs will be saved here
"load_fastICA_results" : False, # If all the FastICA runs already exisit, setting this to True speeds up ICASAR as they don't need to be recomputed.
"load_fastICA_results" : True, # If all the FastICA runs already exisit, setting this to True speeds up ICASAR as they don't need to be recomputed.
"create_all_ifgs_flag" : True, # small signals are hard for ICA to extact from time series, so make it easier by creating all possible long temporal baseline ifgs from the incremental data.
'sica_tica' : 'sica', # controls whether spatial sources or time courses are independent.
'max_n_all_ifgs' : 1000,
Expand All @@ -67,7 +67,7 @@
"ica_param" : (1e-2, 150), # (tolerance, max iterations)
"hdbscan_param" : (100,10), # (min_cluster_size, min_samples) Discussed in more detail in Mcinnes et al. (2017). min_cluster_size sets the smallest collection of points that can be considered a cluster. min_samples sets how conservative the clustering is. With larger values, more points will be considered noise.
"out_folder" : Path('example_spatial_02_outputs_sICA_incremental'), # outputs will be saved here
"load_fastICA_results" : False, # If all the FastICA runs already exisit, setting this to True speeds up ICASAR as they don't need to be recomputed.
"load_fastICA_results" : True, # If all the FastICA runs already exisit, setting this to True speeds up ICASAR as they don't need to be recomputed.
"create_all_ifgs_flag" : False, # small signals are hard for ICA to extact from time series, so make it easier by creating all possible long temporal baseline ifgs from the incremental data.
'max_n_all_ifgs' : 1000,
'sica_tica' : 'sica', # controls whether spatial sources or time courses are independent.
Expand Down
Binary file modified example_spatial_02_outputs_sICA/01_pca_variance_line.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified example_spatial_02_outputs_sICA/02_PCA_sources_correlations.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified example_spatial_02_outputs_sICA/02_PCA_sources_time.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
Binary file added example_spatial_02_outputs_sICA/ICs.kmz
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file modified example_spatial_02_outputs_tICA_incremental/ICASAR_results.pkl
Binary file not shown.
Binary file modified example_spatial_02_outputs_tICA_incremental/ICs.kmz
Binary file not shown.
2 changes: 1 addition & 1 deletion example_temporal_01.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@
"bootstrapping_param" : (200, 0), # (number of runs with bootstrapping, number of runs without bootstrapping) "hdbscan_param" : (35, 10), # (min_cluster_size, min_samples)
"tsne_param" : (30, 12), # (perplexity, early_exaggeration)
"ica_param" : (1e-2, 150), # (tolerance, max iterations)
"hdbscan_param" : (35,10), # (min_cluster_size, min_samples) Discussed in more detail in Mcinnes et al. (2017). min_cluster_size sets the smallest collection of points that can be considered a cluster. min_samples sets how conservative the clustering is. With larger values, more points will be considered noise.
"hdbscan_param" : (50,10), # (min_cluster_size, min_samples) Discussed in more detail in Mcinnes et al. (2017). min_cluster_size sets the smallest collection of points that can be considered a cluster. min_samples sets how conservative the clustering is. With larger values, more points will be considered noise.
"out_folder" : 'example_temporal_01_outputs', # outputs will be saved here
"figures" : "png+window", # if png, saved in a folder as .png. If window, open as interactive matplotlib figures,
"inset_axes_side" : {'x':0.3, 'y':0.1}} # the size of thei inset axes created in the clustering and manifold figure.
Expand Down
Binary file modified example_temporal_01_outputs/01_pca_variance_line.png
Binary file modified example_temporal_01_outputs/02_PCA_sources.png
Binary file modified example_temporal_01_outputs/ICASAR_results.pkl
Binary file not shown.
2 changes: 0 additions & 2 deletions icasar/aux.py
Original file line number Diff line number Diff line change
Expand Up @@ -1172,5 +1172,3 @@ def prepare_legends_for_2d(clusters_by_max_Iq_no_noise, Iq):
'labels' : legend_labels}
return legend_dict


#%%
6 changes: 5 additions & 1 deletion icasar/icasar_funcs.py
Original file line number Diff line number Diff line change
Expand Up @@ -292,12 +292,15 @@ def baselines_from_names(self):
X_mean = ifgs_all.means_space
else:
X_mc = ifgs_dc.mixtures_mc_space # if not, the incremental (daisy chain) of interferograms will be
X_means = ifgs_dc.means_space
X_mean = ifgs_dc.means_space
elif sica_tica == 'tica': # if we're doing temporal ica with spatial data, the mixtures need to be the transpose
X_mc = ifgs_cum.mixtures_mc_time.T # as cumulative and transpose, effectively the time series for each point.
X_mean = ifgs_cum.means_time
else:
X = temporal_data['mixtures_r2']
X_mean = np.mean(X, axis = 1)[:,np.newaxis] # get the mean for each ifg (ie along rows. )
X_mc = X - X_mean # mean centre the data (along rows)




Expand Down Expand Up @@ -443,6 +446,7 @@ def baselines_from_names(self):
inversion_results = bss_components_inversion(S_ica, [X_mc]) # invert to fit the mean centered mixture.
source_residuals = inversion_results[0]['residual'] # how well we fit those
A_ica = inversion_results[0]['tcs'].T # and the time coruses to remake them.
plot_temporal_signals(S_ica, '04_ICASAR_sources', **fig_kwargs)



Expand Down

0 comments on commit ea0e217

Please sign in to comment.