Commit 00da0786 authored by Marco Cristoforetti's avatar Marco Cristoforetti
Browse files

script classifier

parent 5afdcbe6
This diff is collapsed.
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
import sys; sys.path.append('/home/mcristofo/DST/DST')
import os
from DST.config import data_path
import pandas as pd
import numpy as np
import torch
import torch.nn as nn
import time
import math
import torch.utils.data as utils_data
import torch.nn.functional as F
import datetime
from sklearn.metrics import confusion_matrix, matthews_corrcoef
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
np.set_printoptions(suppress=True, precision=3)
torch.manual_seed(21894)
np.random.seed(21894)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(device)
BEFORE = 12
AFTER = 12
dst_data = pd.read_pickle(os.path.join(data_path,'dst.pkl'))
dst_data['ora_round'] = dst_data.ora.apply(lambda x:int(x.split(':')[0]))
dati_agg = dst_data.groupby(['data','ora_round']).agg({
'BX': np.mean,
'BY': np.mean,
'BZ': np.mean,
'FLOW_SPEED': np.mean,
'PROTON_DENSITY': np.mean,
'TEMPERATURE': np.mean,
'PRESSION': np.mean,
'ELETTRIC': np.mean,
'y': np.mean})
dati_agg.reset_index(inplace=True)
dati_agg.sort_values(by = ['data','ora_round'],inplace=True)
dataset = dati_agg.drop(columns = ['data','ora_round']).values
dataset = torch.from_numpy(np.hstack([np.arange(len(dataset)).reshape([-1,1]),dataset]))
last_date_train = dati_agg[dati_agg.data <= datetime.datetime(2008,12,31)].index[-1]
len_valid_test = (len(dataset) - last_date_train)/2
last_date_train/len(dataset), len_valid_test/len(dataset)
data_in = dataset.unfold(0, BEFORE, 1).transpose(2,1)
data_out = dataset[BEFORE:].unfold(0, AFTER, 1).transpose(2,1)
data_in = data_in[:data_out.size(0)]
data_out = data_out[:,:,-1]
data_in.size(), data_out.size()
where_not_nan_in = ~torch.isnan(data_in).any(2, keepdim=True).any(1, keepdim=True).reshape(-1)
data_in = data_in[where_not_nan_in]
data_out = data_out[where_not_nan_in]
where_not_nan_out = ~torch.isnan(data_out).any(1, keepdim=True).reshape(-1)
data_in = data_in[where_not_nan_out]
data_out = data_out[where_not_nan_out]
last_train = np.where(data_in[:,0,0] <= last_date_train)[0][-1] + 1
data_in = data_in[:, :, 1:]
#len_tr = int(data_in.size(0) * 0.6)
n_channels = data_in.size(2)
class MinMaxScaler():
"""
Transform features by scaling each feature to a given range
Features in the last dim
The transformation is given by::
X_std = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0))
X_scaled = X_std * (max - min) + min
where min, max = feature_range.
"""
def __init__(self, feature_range=(0,1)):
self.feature_range = feature_range
def fit(self, X):
X_size = X.size()
X = X.reshape(-1, X_size[-1])
data_min = X.min(axis=0).values
data_max = X.max(axis=0).values
data_range = data_max - data_min
self.scale_ = ((self.feature_range[1] - self.feature_range[0]) / data_range)
self.min_ = self.feature_range[0] - data_min * self.scale_
self.data_min_ = data_min
self.data_max_ = data_max
self.data_range_ = data_range
X = X.reshape(X_size)
return self
def transform(self, X):
X *= self.scale_
X += self.min_
return X
def inverse_transform(self, X):
X -= self.min_
X /= self.scale_
return X
mmScaler = MinMaxScaler((0.1, .9))
mmScaler.fit(data_in[:last_train])
data_in_scaled = data_in.clone()
data_in_scaled = mmScaler.transform(data_in_scaled)
mm_scaler_out = MinMaxScaler((0.1, .9))
mm_scaler_out.fit(data_in[:last_train, :, -1].reshape(-1, data_in.size(1), 1))
data_out_scaled = data_out.clone()
data_out_scaled = mm_scaler_out.transform(data_out_scaled)
dst_levels = [-20,-50,-100]
data_out_c = data_out.clone()
data_out_c[np.where(data_out_c >= dst_levels[0])] = 0
data_out_c[np.where((data_out_c < dst_levels[0]) & (data_out_c >= dst_levels[1]))] = 1
data_out_c[np.where((data_out_c < dst_levels[1]) & (data_out_c >= dst_levels[2]))] = 2
data_out_c[np.where((data_out_c < dst_levels[2]))] = 3
class Dataset(utils_data.Dataset):
def __init__(self, dataset_in, dataset_out, dataset_out_c):
self.dataset_in = dataset_in
self.dataset_out = dataset_out
self.dataset_out_c = dataset_out_c
def __len__(self):
return self.dataset_in.size(0)
def __getitem__(self, idx):
din_src = self.dataset_in[idx]
dout = self.dataset_out[idx]
dout_c = self.dataset_out_c[idx]
return din_src, dout, dout_c
ixs_valid_test = np.arange(int(len_valid_test)) + last_train
np.random.shuffle(ixs_valid_test)
ixs_valid = ixs_valid_test[::2]
ixs_test = ixs_valid_test[1::2]
ixs_tr1 = np.where((mm_scaler_out.inverse_transform(data_out_scaled[:last_train].clone()).min(axis=1)[0]).numpy()<-20)[0]
ixs_tr2 = np.where((mm_scaler_out.inverse_transform(data_out_scaled[:last_train].clone()).min(axis=1)[0]).numpy()>=-20)[0]
BATCH_SIZE=256
np.random.shuffle(ixs_tr2)
ixs_tr2a = ixs_tr2[:10000]
ixs_tr = list(ixs_tr1) + list(ixs_tr2a)
aa = data_out_c[ixs_tr]
weights_c = torch.tensor([len(aa[aa==0])/len(aa[aa==0]), len(aa[aa==0])/len(aa[aa==1]), len(aa[aa==0])/len(aa[aa==2]), len(aa[aa==0])/len(aa[aa==3])]).to(device)#.sqrt()
dataset_tr = Dataset(data_in_scaled[ixs_tr], data_out_scaled[ixs_tr], data_out_c[ixs_tr])
# data_loader_tr = utils_data.DataLoader(dataset_tr, batch_size=BATCH_SIZE, num_workers = 4, shuffle=False, sampler = sampler)
data_loader_tr = utils_data.DataLoader(dataset_tr, batch_size=BATCH_SIZE, num_workers = 4, shuffle=True)
class DSTnet(nn.Module):
def __init__(self, nvars, nhidden_i, nhidden_o, n_out_i, before, after):
super().__init__()
self.nvars = nvars
self.nhidden_i = nhidden_i
self.nhidden_o = nhidden_o
self.before = before
self.after = after
self.n_out_i = n_out_i
self.lstm = nn.LSTM(self.nvars, self.n_out_i, self.nhidden_i, batch_first=True)
self.first_merged_layer = self.n_out_i * self.before
self.bn1 = nn.BatchNorm1d(num_features=self.first_merged_layer)
self.linear_o_1 = nn.Linear(self.first_merged_layer, self.nhidden_o)
self.ln1 = nn.LayerNorm(self.nhidden_o )
self.linear_o_2 = nn.Linear(self.nhidden_o, self.nhidden_o)
self.linear_o_3 = nn.Linear(self.nhidden_o, self.nhidden_o // 2)
self.linear_o_4_c = nn.Linear(self.nhidden_o // 2, self.after * 4)
def init_hidden(self, batch_size):
hidden = torch.randn(self.nhidden_i, batch_size, self.n_out_i).to(device)
cell = torch.randn(self.nhidden_i, batch_size, self.n_out_i).to(device)
return (hidden, cell)
def forward(self, x0):
self.hidden = self.init_hidden(x0.size(0))
x = self.lstm(x0, self.hidden)[0].reshape(x0.shape[0], -1)
x = self.bn1(x)
# x = F.relu(x)
x = F.relu(self.linear_o_1(x))
# x = self.ln1(x)
x = F.dropout(x, 0.2, training=self.training)
x = F.relu(self.linear_o_2(x))
x = F.dropout(x, 0.2, training=self.training)
x = F.relu(self.linear_o_3(x))
x = F.dropout(x, 0.2, training=self.training)
x1 = self.linear_o_4_c(x)
x1 = x1.reshape(x0.size(0) * self.after, 4)
return x1
loss_f = nn.L1Loss()
# loss_fc= nn.CrossEntropyLoss()
loss_fc= nn.CrossEntropyLoss(weight = weights_c)
nhidden_i = 2
nhidden_o = 96
n_out_i = 8
before = BEFORE
nvars = data_in_scaled.shape[-1]
# dst_net = DSTnet(nvars, nhidden_i, nhidden_o, n_out_i, before, AFTER).to(device)
print(dst_net)
num_epochs = 2000
lr = 1e-5
optimizer = torch.optim.Adam(dst_net.parameters(), lr=lr)#, weight_decay=1e-5)
history_tr = np.zeros(num_epochs)
history_valid = np.zeros(num_epochs)
history_ts = np.zeros(num_epochs)
for epoch in range(num_epochs):
start_time = time.time()
for i, batch in enumerate(data_loader_tr):
x = batch[0].float().to(device)
y_r = batch[1].float().to(device)
y_c = batch[2].flatten().long().to(device)
optimizer.zero_grad()
dst_net.train()
out_c = dst_net(x)
loss = loss_fc(out_c, y_c)
loss.backward()
optimizer.step()
dst_net.eval()
out_c = dst_net(data_in_scaled[:last_train].to(device).float())
loss_c_tr = loss_fc(out_c, data_out_c[:last_train].flatten().long().to(device)).item()
out_c = dst_net(data_in_scaled[ixs_valid].to(device).float())
loss_c_val = loss_fc(out_c, data_out_c[ixs_valid].flatten().long().to(device)).item()
out_c = dst_net(data_in_scaled[ixs_test].to(device).float())
loss_c_ts = loss_fc(out_c, data_out_c[ixs_test].flatten().long().to(device)).item()
history_tr[epoch] = loss_c_tr
history_valid[epoch] = loss_c_val
history_ts[epoch] = loss_c_ts
epoch_time = time.time() - start_time
if (epoch % 10 == 0):
print('Epoch %d time = %.2f, tr_c = %.5f, val_c = %.5f, ts_c = %.5f' %
(epoch, epoch_time, loss_c_tr, loss_c_val, loss_c_ts))
if (epoch % 10 == 0):
np.set_printoptions(suppress=True, precision=3)
out_c = dst_net(data_in_scaled[:last_train].to(device).float())
out_cc = F.softmax(out_c, dim=1)
out_cc = out_cc.detach().cpu().numpy()
tp = 11
print(confusion_matrix(data_out_c[:last_train][:,tp].cpu().detach().numpy(), out_cc.argmax(axis=1).reshape(-1,12)[:,tp])/(confusion_matrix(data_out_c[:last_train][:,tp].cpu().detach().numpy(), out_cc.argmax(axis=1).reshape(-1,12)[:,tp]).sum(axis=1)[:, None]))
out_c = dst_net(data_in_scaled[last_train:].to(device).float())
out_cc = F.softmax(out_c, dim=1)
out_cc = out_cc.detach().cpu().numpy()
tp = 11
print(confusion_matrix(data_out_c[last_train:][:,tp].cpu().detach().numpy(), out_cc.argmax(axis=1).reshape(-1,12)[:,tp])/(confusion_matrix(data_out_c[last_train:][:,tp].cpu().detach().numpy(), out_cc.argmax(axis=1).reshape(-1,12)[:,tp]).sum(axis=1)[:, None]))
np.random.shuffle(ixs_tr2)
ixs_tr2a = ixs_tr2[:10000]
ixs_tr = list(ixs_tr1) + list(ixs_tr2a)
dataset_tr = Dataset(data_in_scaled[ixs_tr], data_out_scaled[ixs_tr], data_out_c[ixs_tr])
data_loader_tr = utils_data.DataLoader(dataset_tr, batch_size=BATCH_SIZE, num_workers = 4, shuffle=True)
aa = data_out_c[ixs_tr]
weights_c = torch.tensor([len(aa[aa==0])/len(aa[aa==0]), len(aa[aa==0])/len(aa[aa==1]), len(aa[aa==0])/len(aa[aa==2]), len(aa[aa==0])/len(aa[aa==3])]).to(device)#.sqrt()
loss_fc= nn.CrossEntropyLoss(weight = weights_c)
torch.save(dst_net.state_dict(), os.path.join('/home/mcristofo/DST/models','dst_class.pth'))
np.savetxt(os.path.join('/home/mcristofo/DST/hist','history_tr_class.txt'), history_tr)
np.savetxt(os.path.join('/home/mcristofo/DST/hist','history_valid_class.txt'), history_valid)
np.savetxt(os.path.join('/home/mcristofo/DST/hist','history_ts_class.txt'), history_ts)
import sys; sys.path.append('/home/mcristofo/DST/DST')
import os
from DST.config import data_path
import pandas as pd
import numpy as np
import torch
import torch.nn as nn
import time
import math
import torch.utils.data as utils_data
import torch.nn.functional as F
import datetime
from sklearn.metrics import confusion_matrix, matthews_corrcoef
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
np.set_printoptions(suppress=True, precision=3)
torch.manual_seed(21894)
np.random.seed(21894)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(device)
BEFORE = 12
AFTER = 12
dst_data = pd.read_pickle(os.path.join(data_path,'dst.pkl'))
dst_data['ora_round'] = dst_data.ora.apply(lambda x:int(x.split(':')[0]))
dati_agg = dst_data.groupby(['data','ora_round']).agg({
'BX': np.mean,
'BY': np.mean,
'BZ': np.mean,
'FLOW_SPEED': np.mean,
'PROTON_DENSITY': np.mean,
'TEMPERATURE': np.mean,
'PRESSION': np.mean,
'ELETTRIC': np.mean,
'y': np.mean})
dati_agg.reset_index(inplace=True)
dati_agg.sort_values(by = ['data','ora_round'],inplace=True)
dataset = dati_agg.drop(columns = ['data','ora_round']).values
dataset = torch.from_numpy(np.hstack([np.arange(len(dataset)).reshape([-1,1]),dataset]))
last_date_train = dati_agg[dati_agg.data <= datetime.datetime(2008,12,31)].index[-1]
len_valid_test = (len(dataset) - last_date_train)/2
last_date_train/len(dataset), len_valid_test/len(dataset)
data_in = dataset.unfold(0, BEFORE, 1).transpose(2,1)
data_out = dataset[BEFORE:].unfold(0, AFTER, 1).transpose(2,1)
data_in = data_in[:data_out.size(0)]
data_out = data_out[:,:,-1]
data_in.size(), data_out.size()
where_not_nan_in = ~torch.isnan(data_in).any(2, keepdim=True).any(1, keepdim=True).reshape(-1)
data_in = data_in[where_not_nan_in]
data_out = data_out[where_not_nan_in]
where_not_nan_out = ~torch.isnan(data_out).any(1, keepdim=True).reshape(-1)
data_in = data_in[where_not_nan_out]
data_out = data_out[where_not_nan_out]
last_train = np.where(data_in[:,0,0] <= last_date_train)[0][-1] + 1
data_in = data_in[:, :, 1:]
#len_tr = int(data_in.size(0) * 0.6)
n_channels = data_in.size(2)
class MinMaxScaler():
"""
Transform features by scaling each feature to a given range
Features in the last dim
The transformation is given by::
X_std = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0))
X_scaled = X_std * (max - min) + min
where min, max = feature_range.
"""
def __init__(self, feature_range=(0,1)):
self.feature_range = feature_range
def fit(self, X):
X_size = X.size()
X = X.reshape(-1, X_size[-1])
data_min = X.min(axis=0).values
data_max = X.max(axis=0).values
data_range = data_max - data_min
self.scale_ = ((self.feature_range[1] - self.feature_range[0]) / data_range)
self.min_ = self.feature_range[0] - data_min * self.scale_
self.data_min_ = data_min
self.data_max_ = data_max
self.data_range_ = data_range
X = X.reshape(X_size)
return self
def transform(self, X):
X *= self.scale_
X += self.min_
return X
def inverse_transform(self, X):
X -= self.min_
X /= self.scale_
return X
mmScaler = MinMaxScaler((0.1, .9))
mmScaler.fit(data_in[:last_train])
data_in_scaled = data_in.clone()
data_in_scaled = mmScaler.transform(data_in_scaled)
mm_scaler_out = MinMaxScaler((0.1, .9))
mm_scaler_out.fit(data_in[:last_train, :, -1].reshape(-1, data_in.size(1), 1))
data_out_scaled = data_out.clone()
data_out_scaled = mm_scaler_out.transform(data_out_scaled)
dst_levels = [-20,-50,-100]
data_out_c = data_out.clone()
data_out_c[np.where(data_out_c >= dst_levels[0])] = 0
data_out_c[np.where((data_out_c < dst_levels[0]) & (data_out_c >= dst_levels[1]))] = 1
data_out_c[np.where((data_out_c < dst_levels[1]) & (data_out_c >= dst_levels[2]))] = 2
data_out_c[np.where((data_out_c < dst_levels[2]))] = 3
class Dataset(utils_data.Dataset):
def __init__(self, dataset_in, dataset_out, dataset_out_c):
self.dataset_in = dataset_in
self.dataset_out = dataset_out
self.dataset_out_c = dataset_out_c
def __len__(self):
return self.dataset_in.size(0)
def __getitem__(self, idx):
din_src = self.dataset_in[idx]
dout = self.dataset_out[idx]
dout_c = self.dataset_out_c[idx]
return din_src, dout, dout_c
ixs_valid_test = np.arange(int(len_valid_test)) + last_train
np.random.shuffle(ixs_valid_test)
ixs_valid = ixs_valid_test[::2]
ixs_test = ixs_valid_test[1::2]
ixs_tr1 = np.where((mm_scaler_out.inverse_transform(data_out_scaled[:last_train].clone()).min(axis=1)[0]).numpy()<-20)[0]
ixs_tr2 = np.where((mm_scaler_out.inverse_transform(data_out_scaled[:last_train].clone()).min(axis=1)[0]).numpy()>=-20)[0]
BATCH_SIZE=256
np.random.shuffle(ixs_tr2)
ixs_tr2a = ixs_tr2[:10000]
ixs_tr = list(ixs_tr1) + list(ixs_tr2a)
aa = data_out_c[ixs_tr]
weights_c = torch.tensor([len(aa[aa==0])/len(aa[aa==0]), len(aa[aa==0])/len(aa[aa==1]), len(aa[aa==0])/len(aa[aa==2]), len(aa[aa==0])/len(aa[aa==3])]).to(device)#.sqrt()
dataset_tr = Dataset(data_in_scaled[ixs_tr], data_out_scaled[ixs_tr], data_out_c[ixs_tr])
# data_loader_tr = utils_data.DataLoader(dataset_tr, batch_size=BATCH_SIZE, num_workers = 4, shuffle=False, sampler = sampler)
data_loader_tr = utils_data.DataLoader(dataset_tr, batch_size=BATCH_SIZE, num_workers = 4, shuffle=True)
class DSTnet(nn.Module):
def __init__(self, nvars, nhidden_i, nhidden_o, n_out_i, before, after):
super().__init__()
self.nvars = nvars
self.nhidden_i = nhidden_i
self.nhidden_o = nhidden_o
self.before = before
self.after = after
self.n_out_i = n_out_i
self.lstm = nn.LSTM(self.nvars, self.n_out_i, self.nhidden_i, batch_first=True)
self.first_merged_layer = self.n_out_i * self.before
self.bn1 = nn.BatchNorm1d(num_features=self.first_merged_layer)
self.linear_o_1 = nn.Linear(self.first_merged_layer, self.nhidden_o)
self.ln1 = nn.LayerNorm(self.nhidden_o )
self.linear_o_2 = nn.Linear(self.nhidden_o, self.nhidden_o)
self.linear_o_3 = nn.Linear(self.nhidden_o, self.nhidden_o // 2)
self.linear_o_4_c = nn.Linear(self.nhidden_o // 2, self.after * 4)
def init_hidden(self, batch_size):
hidden = torch.randn(self.nhidden_i, batch_size, self.n_out_i).to(device)
cell = torch.randn(self.nhidden_i, batch_size, self.n_out_i).to(device)
return (hidden, cell)
def forward(self, x0):
self.hidden = self.init_hidden(x0.size(0))
x = self.lstm(x0, self.hidden)[0].reshape(x0.shape[0], -1)
x = self.bn1(x)
# x = F.relu(x)
x = F.relu(self.linear_o_1(x))
# x = self.ln1(x)
x = F.dropout(x, 0.2, training=self.training)
x = F.relu(self.linear_o_2(x))
x = F.dropout(x, 0.2, training=self.training)
x = F.relu(self.linear_o_3(x))
x = F.dropout(x, 0.2, training=self.training)
x1 = self.linear_o_4_c(x)
x1 = x1.reshape(x0.size(0) * self.after, 4)
return x1
loss_f = nn.L1Loss()
# loss_fc= nn.CrossEntropyLoss()
loss_fc= nn.CrossEntropyLoss(weight = weights_c)
nhidden_i = 2
nhidden_o = 96
n_out_i = 8
before = BEFORE
nvars = data_in_scaled.shape[-1]
# dst_net = DSTnet(nvars, nhidden_i, nhidden_o, n_out_i, before, AFTER).to(device)
print(dst_net)
num_epochs = 2000
lr = 1e-5
optimizer = torch.optim.Adam(dst_net.parameters(), lr=lr)#, weight_decay=1e-5)
history_tr = np.zeros(num_epochs)
history_valid = np.zeros(num_epochs)
history_ts = np.zeros(num_epochs)
for epoch in range(num_epochs):
start_time = time.time()
for i, batch in enumerate(data_loader_tr):
x = batch[0].float().to(device)
y_r = batch[1].float().to(device)
y_c = batch[2].flatten().long().to(device)
optimizer.zero_grad()
dst_net.train()