数据
df.to_csv('traffic1.txt', header=None, index=None)
df
数据为每小时记录一次的交通流量数据,每周有几天出现早高峰。
dataset
import torch
import numpy as np;
from torch.autograd import Variable
def normal_std(x):
return x.std() * np.sqrt((len(x) - 1.)/(len(x)))
class Data_utility(object):
# train and valid is the ratio of training set and validation set. test = 1 - train - valid
def __init__(self, file_name, train, valid, cuda, horizon, window, normalize = 2):
self.cuda = cuda;
self.P = window;
self.h = horizon
fin = open(file_name);
self.rawdat = np.loadtxt(fin,delimiter=',')
if(len(self.rawdat.shape)) == 1:
self.rawdat = self.rawdat.reshape(len(self.rawdat), -1)
self.dat = np.zeros(self.rawdat.shape)
self.n, self.m = self.dat.shape;
self.normalize = 2
self.scale = np.ones(self.m);
self._normalized(normalize);
self._split(int(train * self.n), int((train+valid) * self.n), self.n);
self.scale = torch.from_numpy(self.scale).float();
tmp = self.test[1] * self.scale.expand(self.test[1].size(0), self.m);
if self.cuda:
self.scale = self.scale.cuda();
self.scale = Variable(self.scale);
self.rse = normal_std(tmp);
self.rae = torch.mean(torch.abs(tmp - torch.mean(tmp)));
def _normalized(self, normalize):
#normalized by the maximum value of entire matrix.
if (normalize == 0):
self.dat = self.rawdat
if (normalize == 1):
self.dat = self.rawdat / np.max(self.rawdat);
#normlized by the maximum value of each row(sensor).
if (normalize == 2):
for i in range(self.m):
self.scale[i] = np.max(np.abs(self.rawdat[:,i]));
self.dat[:,i] = self.rawdat[:,i] / np.max(np.abs(self.rawdat[:,i]));
def _split(self, train, valid, test):
train_set = range(self.P+self.h-1, train);
valid_set = range(train, valid);
test_set = range(valid, self.n);
self.train = self._batchify(train_set, self.h);
self.valid = self._batchify(valid_set, self.h);
self.test = self._batchify(test_set, self.h);
def _batchify(self, idx_set, horizon):
n = len(idx_set);
X = torch.zeros((n,self.P,self.m));
Y = torch.zeros((n,self.m));
for i in range(n):
end = idx_set[i] - self.h + 1;
start = end - self.P;
X[i,:,:] = torch.from_numpy(self.dat[start:end, :]);
Y[i,:] = torch.from_numpy(self.dat[idx_set[i], :]);
return [X, Y];
def get_batches(self, inputs, targets, batch_size, shuffle=True):
length = len(inputs)
if shuffle:
index = torch.randperm(length)
else:
index = torch.LongTensor(range(length))
start_idx = 0
while (start_idx < length):
end_idx = min(length, start_idx + batch_size)
excerpt = index[start_idx:end_idx]
X, Y= inputs[excerpt], targets[excerpt]
if (self.cuda):
X, Y = X.cuda(), Y.cuda()
yield Variable(X), Variable(Y)
start_idx += batch_size
data = Data_utility(file_name='traffic1.txt', train=0.6, valid=0.2, cuda=False, horizon=12, window=24*7, normalize = 2)
print(data.train[0].shape, data.train[1].shape) # torch.Size([10347, 168, 1]) torch.Size([10347, 1])
window = data.train[0].shape[1]
n_val = data.train[0].shape[2]
optimizer
import math
import torch.optim as optim
class Optim(object):
def _makeOptimizer(self):
if self.method == 'sgd':
self.optimizer = optim.SGD(self.params, lr=self.lr)
elif self.method == 'adagrad':
self.optimizer = optim.Adagrad(self.params, lr=self.lr)
elif self.method == 'adadelta':
self.optimizer = optim.Adadelta(self.params, lr=self.lr)
elif self.method == 'adam':
self.optimizer = optim.Adam(self.params, lr=self.lr)
else:
raise RuntimeError("Invalid optim method: " + self.method)
def __init__(self, params, method, lr, max_grad_norm, lr_decay=1, start_decay_at=None):
self.params = list(params) # careful: params may be a generator
self.last_ppl = None
self.lr = lr
self.max_grad_norm = max_grad_norm
self.method = method
self.lr_decay = lr_decay
self.start_decay_at = start_decay_at
self.start_decay = False
self._makeOptimizer()
def step(self):
# Compute gradients norm.
grad_norm = 0
for param in self.params:
grad_norm += math.pow(param.grad.data.norm(), 2)
grad_norm = math.sqrt(grad_norm)
if grad_norm > 0:
shrinkage = self.max_grad_norm / grad_norm
else:
shrinkage = 1.
for param in self.params:
if shrinkage < 1:
param.grad.data.mul_(shrinkage)
self.optimizer.step()
return grad_norm
# decay learning rate if validation performance does not improve or we hit the start_decay_at limit
def updateLearningRate(self, ppl, epoch):
if self.start_decay_at is None or epoch <= self.start_decay_at:
return
if self.last_ppl is not None and ppl > self.last_ppl:
self.start_decay = True
if self.start_decay:
self.lr = self.lr * self.lr_decay
# print("Decaying learning rate to %g" % self.lr)
#only decay for one epoch
self.start_decay = False
self.last_ppl = ppl
self._makeOptimizer()
model
import torch
import torch.nn as nn
import torch.nn.functional as F
class Model(nn.Module):
def __init__(self, n_val, window, hidRNN):
super(Model, self).__init__()
self.use_cuda = True
self.P = window # 输入窗口大小
self.m = n_val # 列数,变量数
self.hidR = hidRNN
self.GRU = nn.GRU(self.m, self.hidR)
self.linear = nn.Linear(self.hidR, self.m)
def forward(self, x):
# x: [batch, window, n_val]
# batch_size = x.shape[0]
# x_flat = x.view(batch_size, -1)
x1 = x.permute(1, 0, 2).contiguous() # x1: [window, batch, n_val]
_, h = self.GRU(x1) # r: [1, batch, hidRNN]
h = torch.squeeze(h,0) # r: [batch, hidRNN]
res = self.linear(h) # res: [batch, n_val]
return res
train
def train(data, X, Y, model, criterion, optim, batch_size):
model.train()
total_loss = 0
n_samples = 0
for X, Y in data.get_batches(X, Y, batch_size, True):
model.zero_grad();
output = model(X);
scale = data.scale.expand(output.size(0), data.m)
loss = criterion(output * scale, Y * scale);
loss.backward();
grad_norm = optim.step();
total_loss += loss.item();
n_samples += (output.size(0) * data.m);
return total_loss / n_samples
def evaluate(data, X, Y, model, evaluateL2, evaluateL1, batch_size):
model.eval()
total_loss = 0
total_loss_l1 = 0
n_samples = 0
predict = None
test = None
for X, Y in data.get_batches(X, Y, batch_size, False):
output = model(X)
if predict is None:
predict = output
test = Y
else:
predict = torch.cat((predict,output))
test = torch.cat((test, Y))
scale = data.scale.expand(output.size(0), data.m)
total_loss += evaluateL2(output * scale, Y * scale).item()
total_loss_l1 += evaluateL1(output * scale, Y * scale).item()
n_samples += (output.size(0) * data.m)
rse = math.sqrt(total_loss / n_samples)/data.rse
rae = (total_loss_l1/n_samples)/data.rae
predict = predict.data.cpu().numpy();
Ytest = test.data.cpu().numpy();
sigma_p = (predict).std(axis = 0);
sigma_g = (Ytest).std(axis = 0);
mean_p = predict.mean(axis = 0)
mean_g = Ytest.mean(axis = 0)
index = (sigma_g!=0);
correlation = ((predict - mean_p) * (Ytest - mean_g)).mean(axis = 0)/(sigma_p * sigma_g);
correlation = (correlation[index]).mean();
return rse, rae, correlation, predict
model = Model(n_val, window, 32);
nParams = sum([p.nelement() for p in model.parameters()])
print('* number of parameters: %d' % nParams)
criterion = nn.MSELoss(reduction='sum')
evaluateL2 = nn.MSELoss(reduction='sum')
evaluateL1 = nn.L1Loss(reduction='sum')
criterion = criterion.cuda()
evaluateL1 = evaluateL1.cuda()
evaluateL2 = evaluateL2.cuda()
optimizer = Optim(
model.parameters(), 'adam', lr=0.01, max_grad_norm=10, start_decay_at = 10, lr_decay = 0.9
)
batch_size=128
epochs =10
best_val = 10000000
save = 'model.pt'
print('begin training')
import time
for epoch in range(1, epochs):
epoch_start_time = time.time()
train_loss = train(data, data.train[0], data.train[1], model, criterion, optimizer, batch_size)
val_loss, val_rae, val_corr, _ = evaluate(data, data.valid[0], data.valid[1], model, evaluateL2, evaluateL1, batch_size);
print('| end of epoch {:3d} | time: {:5.2f}s | train_loss {:5.4f} | valid rse {:5.4f} | valid rae {:5.4f} | valid corr {:5.4f} | lr {:5.4f}'
.format(epoch, (time.time() - epoch_start_time), train_loss, val_loss, val_rae, val_corr, optimizer.lr))
# Save the model if the validation loss is the best we've seen so far.
if val_loss < best_val:
with open(save, 'wb') as f:
torch.save(model, f)
best_val = val_loss
if epoch % 5 == 0:
test_acc, test_rae, test_corr, _ = evaluate(data, data.test[0], data.test[1], model, evaluateL2, evaluateL1, batch_size);
print ("test rse {:5.4f} | test rae {:5.4f} | test corr {:5.4f}".format(test_acc, test_rae, test_corr))
optimizer.updateLearningRate(val_loss, epoch)
* number of parameters: 3393
begin training
| end of epoch 1 | time: 21.82s | train_loss 0.0036 | valid rse 0.6420 | valid rae 0.7820 | valid corr 0.4746 | lr 0.0100
| end of epoch 2 | time: 22.95s | train_loss 0.0024 | valid rse 0.5568 | valid rae 0.5982 | valid corr 0.5768 | lr 0.0100
| end of epoch 3 | time: 25.06s | train_loss 0.0021 | valid rse 0.5137 | valid rae 0.5141 | valid corr 0.6616 | lr 0.0100
| end of epoch 4 | time: 23.82s | train_loss 0.0020 | valid rse 0.4988 | valid rae 0.4964 | valid corr 0.7043 | lr 0.0100
| end of epoch 5 | time: 23.95s | train_loss 0.0021 | valid rse 0.5104 | valid rae 0.4532 | valid corr 0.6580 | lr 0.0100
test rse 0.7352 | test rae 0.6122 | test corr 0.6882
| end of epoch 6 | time: 23.32s | train_loss 0.0020 | valid rse 0.4880 | valid rae 0.4146 | valid corr 0.7130 | lr 0.0100
| end of epoch 7 | time: 23.69s | train_loss 0.0020 | valid rse 0.4840 | valid rae 0.4461 | valid corr 0.6972 | lr 0.0100
| end of epoch 8 | time: 24.21s | train_loss 0.0019 | valid rse 0.4765 | valid rae 0.4374 | valid corr 0.7112 | lr 0.0100
| end of epoch 9 | time: 23.89s | train_loss 0.0019 | valid rse 0.4779 | valid rae 0.4434 | valid corr 0.7292 | lr 0.0100
test_acc, test_rae, test_corr, pred = evaluate(data, data.test[0], data.test[1], model, evaluateL2, evaluateL1, batch_size);
print ("test rse {:5.4f} | test rae {:5.4f} | test corr {:5.4f}".format(test_acc, test_rae, test_corr))
test rse 0.6917 | test rae 0.5760 | test corr 0.7227
truth = data.test[1].numpy()
for i in range(truth.shape[1]):
plt.figure(figsize=(10,5))
plt.plot(truth[:24*7*4,i], label='ground truth')
plt.plot(pred[:24*7*4,i], label='prediction')
plt.legend()
plt.show()