-
Notifications
You must be signed in to change notification settings - Fork 1
/
run_param_sweep_LASSO.py
124 lines (108 loc) · 4.97 KB
/
run_param_sweep_LASSO.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
from synthmix import SynthMix
from endmember_class import spec_lib
from unmix import LASSO_unmix
import numpy as np
import pandas as pd
import copy
import matplotlib.pyplot as plt
import random
import os
import datetime
now = datetime.datetime.now()
# create directory for today's date
today = now.strftime("%Y-%m-%d")
todays_results = "results/" + today
try:
os.mkdir(todays_results)
except OSError:
print("Creation of the directory %s failed" % todays_results)
else:
print("Successfully created the directory %s " % todays_results)
result_path = todays_results + "/"
# experiment parameters
lambdas = np.logspace(-7, -1, 70) # weights of L1 regularization term
K = 100 # number of mixtures to test
gauss_dev = 2.5e-4 # standard deviation of gaussian noise added to mixture
min_endmemb = 0.01 # minimum endmember percentage
thresh = 0.01 # threshold for determining presence of endmembers
# initialize metrics
Error_L1_mean = [] # average L1 error on abundance prediction
Error_L2_mean = [] # average L2 error on abundance prediction
RMS_noisy_mean = [] # RMS error on signal reconstruction - noisy signal averaged over K mixtures
RMS_true_mean = [] # RMS error on signal reconstruction - true signal averaged over K mixtures
accuracy_mean = []
precision_mean = []
recall_mean = []
SNR_mean = [] # SNR averaged over K mixtures
results = pd.DataFrame()
# sweep over the parameter lambda
for lam in lambdas:
RMS_noisy_total = 0
RMS_true_total = 0
Error_L1_total = 0
Error_L2_total = 0
SNR_total = 0
true_positive_total = 0
true_negative_total = 0
accuracy_total = 0
precision_total = 0
recall_total = 0
# initialize experiment
np.random.seed(1)
random.seed(1)
endmembs = spec_lib("asu",
ascii_spectra="input/endmemb_libs/rogers_tes73.txt",
meta_csv="input/endmemb_libs/rogers_meta.csv")
# endmembs = endmembs.at_most(1)
synth = SynthMix(thresh=min_endmemb)
metrics = []
N = endmembs.spectra.shape[0] # number of spectral channels
M = endmembs.spectra.shape[1] # number of endmembers
# iterate over K mixtures
for i in range(K):
# create synthetic mixture
mixture = synth.mars_mix(endmembs)
mixture_noisy = mixture.perturb(method="gauss", deviation=gauss_dev)
# predict abundances
x = LASSO_unmix(endmembs.spectra, mixture_noisy.spectra, lam=lam)
x = x / sum(x)
x_presence = np.zeros(x.shape)
x_presence[np.nonzero(x > thresh)] = 1
# calculate loss metrics
idx_pos_truth = np.nonzero(mixture.presence > 0)[0]
idx_neg_truth = np.nonzero(mixture.presence == 0)[0]
RMS_noisy_total += np.asscalar((1.0 / N) * sum((mixture_noisy.spectra - np.matmul(endmembs.spectra, x)) ** 2))
RMS_true_total += np.asscalar((1.0 / N) * sum((mixture.spectra - np.matmul(endmembs.spectra, x)) ** 2))
Error_L1_total += np.asscalar(sum(abs(x - mixture.proportions.transpose())))
Error_L2_total += np.asscalar(sum((x - mixture.proportions.transpose()) ** 2))
accuracy_total += sum(mixture.presence == x_presence) / len(mixture.presence)
precision_total += sum(mixture.presence[idx_pos_truth] == x_presence[idx_pos_truth]) / sum(x_presence)
recall_total += sum(mixture.presence[idx_pos_truth] == x_presence[idx_pos_truth]) / sum(mixture.presence)
SNR_total += np.asscalar(sum((mixture.spectra - mixture.spectra.max())**2) / sum((mixture.spectra - mixture_noisy.spectra)**2))
# store average metrics
Error_L1_mean.append((1.0/K) * Error_L1_total)
Error_L2_mean.append((1.0/K) * Error_L2_total)
accuracy_mean.append((1.0/K) * accuracy_total)
precision_mean.append((1.0/K) * precision_total)
recall_mean.append((1.0/K) * recall_total)
RMS_noisy_mean.append((1.0/K) * RMS_noisy_total)
RMS_true_mean.append((1.0/K) * RMS_true_total)
SNR_mean.append((1.0/K) * SNR_total)
print("lambda: %f, Error_L1_mean: %f" % (lam, Error_L1_mean[-1]))
# store results
results = pd.DataFrame({"lambda": lambdas,
"thresh": thresh,
"K": K,
"SNR": SNR_mean,
"Error_L1_mean": Error_L1_mean,
"Error_L2_mean": Error_L2_mean,
"accuracy_mean": accuracy_mean,
"precision_mean": precision_mean,
"recall_mean": recall_mean,
"RMS_noisy_mean": RMS_noisy_mean,
"RMS_true_mean": RMS_true_mean,
"Deviation": gauss_dev})
i = 0
while os.path.exists(result_path + today + "_LASSO_sweep_results%s.csv" % i):
i += 1
results.to_csv(result_path + today + "_LASSO_sweep_results%s.csv" % i)