这个稀疏自编码器就是在编码器的基础上，加上稀疏约束。
所谓稀疏约束就是尽量使隐藏层中激活的神经元是稀少的；
代码里我们算出隐藏层中神经元激活值的平均值q = torch.nn.functional.softmax(encoder_out, dim=1)；
p是我们希望的平均激活值，因为是希望神经元稀疏，所以代码里设置了一个比较小的值：0.2；
再算出q与p的分布差，又称为kl散度。 _kl = criterion_(p, q)；
给损失函数加上这个kl散度： loss += _beta * _kl，这样的目的是使p尽量接近于q，也就是使隐藏层中神经元激活值的平均值接近于0.2，以此间接达到了使神经元稀疏的作用。
以下给出代码，分为训练部分和探测部分：
1、训练部分
'''
训练系数自动编码器用于高光谱异常探测
'''
import torch
from torch import nn, optim
from torch.autograd import Variable
from torchvision import transforms
from torch.utils.data import DataLoader
from torchvision.utils import save_image
import os
import numpy as np
# 加载数据集
def get_data():
'''
函数作用：读取数据集
return：dataload与img（高光谱图像）
'''
from scipy import io
#读取数据
dataset_name='sandiego_plane.mat'
load=io.loadmat(dataset_name)
img = load['data']
gt = load['map']
height,width,spectrum=img.shape
#这一段是预处理数据
img_data=img.reshape((height*width,spectrum))#改变数据形状，使其适合送入dataset
gt_data=gt.reshape((height*width))
img_data=img_data.astype(np.float32)#数据类型转换
# 将像素点转换到[-1, 1]之间，使得输入变成一个比较对称的分布，训练容易收敛
img_data=(img_data-(0.5*img_data.max()))/ (0.5*img_data.max())
#设置dataset与dataload
from datasets import Dataset #这个datasets是自己写的代码文件，见第三段代码
from torch.utils import data
train_dataset = Dataset(img_data, gt_data,dataset_name=dataset_name)
train_loader = DataLoader(train_dataset, shuffle=True, batch_size=batch_size, drop_last=True)
return train_loader,img
def to_img(x):
x = (x + 1.) * 0.5
x = x.clamp(0, 1)
x = x.view(x.size(0), 1, 28, 28)
return x
class autoencoder(nn.Module):
def __init__(self,in_dim,hidden_size,):
super(autoencoder, self).__init__()
self.encoder = nn.Sequential(nn.Linear(in_dim, 128),
nn.ReLU(True),
nn.Linear(128, 64),
nn.ReLU(True),
nn.Linear(64, hidden_size),
nn.ReLU(True)
)
self.decoder = nn.Sequential(nn.Linear(hidden_size, 64),
nn.ReLU(True),
nn.Linear(64, 128),
nn.ReLU(True),
nn.Linear(128, in_dim),
nn.Tanh())
def forward(self, x):
encode = self.encoder(x)
decode = self.decoder(encode)
return encode, decode
def KL_devergence(p, q):
"""
Calculate the KL-divergence of (p,q)，计算p和q之间的KL散度
:param p:期望激活值扩充为hidden_size的tensor
:param q:q为隐藏层结点激活后的输出值
:return:
"""
q = torch.nn.functional.softmax(q, dim=0)# 首先要用softmax对隐藏层结点输出进行归一化处理
q = torch.sum(q, dim=0)/batch_size # dim:缩减的维度,q的第一维是batch维,即大小为batch_size大小,此处是将第j个神经元在batch_size个输入下所有的输出取平均
s1 = torch.sum(p*torch.log(p/q))
s2 = torch.sum((1-p)*torch.log((1-p)/(1-q)))
return s1+s2
if __name__ == "__main__":
# 超参数设置
batch_size = 100
lr = 1e-4
weight_decay = 1e-5
epoches = 100
expect_tho = 0.2
# 读取数据集
train_data,img = get_data()
height, width, spectrum = img.shape
#建立模型
hidden_size=12#隐藏层维度
model = autoencoder(in_dim=spectrum,hidden_size=hidden_size)#参数分别为输入维度与隐藏层维度
criterion = nn.MSELoss()
optimizer = optim.Adam(model.parameters(), lr=lr)#, weight_decay=weight_decay
scheduler=torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode='min', factor=0.3, patience=4, verbose=True,
threshold=0.01, threshold_mode='rel', min_lr=0,
eps=1e-15)
if torch.cuda.is_available():
model.cuda()
# 定义期望平均激活值和KL散度的权重
tho_tensor = torch.FloatTensor([expect_tho for _ in range(hidden_size)])
if torch.cuda.is_available():
tho_tensor = tho_tensor.cuda()
_beta = 1 #通过β来控制KL散度所占的比重.
for epoch in range(epoches):
loss_of_epoch=0
for img, _ in train_data:
# img = Variable(img.cuda())
# forward
encoder_out, output = model(img)
loss = criterion(output, img)
# import torch.functional as F
import torch.nn as F
p = torch.nn.functional.log_softmax(tho_tensor, dim=0)
# q = torch.nn.functional.softmax(encoder_out, dim=1)
q = torch.nn.functional.softmax(encoder_out, dim=1)
# x = F.LogSoftmax(tho_tensor)
# y = F.Softmax(encoder_out, dim=1)
criterion_ = nn.KLDivLoss(reduction = 'batchmean')
_kl = criterion_(p, q)
'''from scipy import stats
# torch.nn.log_softmax
p=tho_tensor.log_softmax().detach().numpy()
q=encoder_out.softmax(y, dim=1).detach().numpy()+1e-10
_kl = stats.entropy(p,q)#计算kl散度
_kl = np.sum(_kl)
_kl = torch.tensor(_kl)'''
# _kl = KL_devergence(tho_tensor, encoder_out)#计算kl散度
loss += _beta * _kl #给损失函数加上kl散度，完成稀疏性约束
# backward
optimizer.zero_grad()
loss.backward()
optimizer.step()
loss_of_epoch = loss_of_epoch+loss.item() # 防止梯度回传
loss_of_epoch/=batch_size#计算一个epoch的平均损失
scheduler.step(loss_of_epoch)
print("epoch: {}, loss is {}".format((epoch+1), loss.data))
for param_group in optimizer.param_groups:
print(param_group['lr'])
torch.save(model.state_dict(), './autoencoder_state.pth')
model(img)2、探测部分
import torch
import numpy as np
from torch import nn
from train_SAE_pytorch import autoencoder,get_data #这里引用的就是上一个代码中的函数
from einops import reduce
def make_AUC(x_test,label):
import sklearn
from sklearn import metrics
x_test=(x_test-np.min(x_test))/(np.max(x_test)-np.min(x_test))
fpr, tpr, _ = sklearn.metrics.roc_curve(label.ravel(), x_test.ravel())
auc = metrics.roc_auc_score(label.ravel(), x_test.ravel(), average='macro')
'''
绘ROC图
'''
import matplotlib.pyplot as plt
plt.figure(figsize=(8, 7), dpi=80, facecolor='w') # dpi:每英寸长度的像素点数；facecolor 背景颜色
plt.xlim((-0.01, 1.02)) # x,y 轴刻度的范围
plt.ylim((-0.01, 1.02))
plt.xticks(np.arange(0, 1.1, 0.1)) # 绘制刻度
plt.yticks(np.arange(0, 1.1, 0.1))
plt.plot(fpr, tpr, 'r-', lw=2, label='AUC=%.4f' % auc) # 绘制AUC 曲线
plt.legend(loc='lower right') # 设置显示标签的位置
plt.xlabel('False Positive Rate', fontsize=14) # 绘制x,y 坐标轴对应的标签
plt.ylabel('True Positive Rate', fontsize=14)
plt.grid(b=True, ls=':') # 绘制网格作为底板;b是否显示网格线；ls表示line style
plt.title(u'DecisionTree ROC curve And AUC', fontsize=18) # 打印标题
plt.show()
if __name__ == "__main__":
import os
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
#加载数据
from scipy import io
dataset_name = 'sandiego_plane.mat'
load = io.loadmat(dataset_name)
img = load['data']
gt = load['map']
height, width, spectrum = img.shape
#建立模型
hidden_size=12#隐藏层维度
#创建模型
model = autoencoder(in_dim=spectrum,hidden_size=hidden_size)
#加载上次训练的模型
try:
model_path='autoencoder_state.pth'
model.load_state_dict(torch.load(model_path,map_location=torch.device('cpu')))
except FileNotFoundError as e:
print(e)
print(f'Loading of model succesful.')
# 这一段是预处理数据
img_data = img.reshape((height * width, spectrum)) # 改变数据形状，使其适合送入dataset
gt_data = gt.reshape((height * width))
img_data = img_data.astype(np.float32) # 数据类型转换
# 将像素点转换到[-1, 1]之间，使得输入变成一个比较对称的分布，训练容易收敛
img_data = (img_data - (0.5 * img_data.max())) / (0.5 * img_data.max())
result=model(torch.tensor(img_data))[1]#转换为tensor才能传入神经网络中 #result为重建结果
#计算异常图， 使用detach剥离梯度属性然后才能使用numpy（）转换数据类型从而进行计算
Anomaly_map=(result.detach().numpy()-img_data)**2 #这里是计算重建误差
Anomaly_map=reduce(Anomaly_map, 'i vec -> i', 'mean')
Anomaly_map=Anomaly_map.reshape((height , width))
#显示
import matplotlib.pyplot as plt
plt.imshow(Anomaly_map, cmap='gray')
plt.show()
make_AUC(Anomaly_map, gt3、用到的函数，命名为 datasets.py 以供调用
import numpy as np
import torch
import torch.utils
import torch.utils.data
import os
from tqdm import tqdm
from sklearn import preprocessing
try:
# Python 3
from urllib.request import urlretrieve
except ImportError:
# Python 2
from urllib import urlretrieve
class TqdmUpTo(tqdm):
"""Provides `update_to(n)` which uses `tqdm.update(delta_n)`."""
def update_to(self, b=1, bsize=1, tsize=None):
"""
b : int, optional
Number of blocks transferred so far [default: 1].
bsize : int, optional
Size of each block (in tqdm units) [default: 1].
tsize : int, optional
Total size (in tqdm units). If [default: None] remains unchanged.
"""
if tsize is not None:
self.total = tsize
self.update(b * bsize - self.n) # will also set self.n = b * bsize
class Dataset(torch.utils.data.Dataset):
""" Generic class for a hyperspectral scene """
def __init__(self, data, gt,dataset_name, **hyperparams):
"""
Args:
data: 3D hyperspectral image
gt: 2D array of labels
patch_size: int, size of the spatial neighbourhood
center_pixel: bool, set to True to consider only the label of the
center pixel
data_augmentation: bool, set to True to perform random flips
supervision: 'full' or 'semi' supervised algorithms
"""
super(Dataset, self).__init__()
self.data = data
self.label = gt
#将label中大于非零值全部置为1，等同于将gt中所有非零值变成1
mask = np.ones_like(gt)
mask[gt == 0] = 0
self.label =mask
@staticmethod
def flip(*arrays):
horizontal = np.random.random() > 0.5
vertical = np.random.random() > 0.5
if horizontal:
arrays = [np.fliplr(arr) for arr in arrays]
if vertical:
arrays = [np.flipud(arr) for arr in arrays]
return arrays
def __len__(self): # 这里返回的是数据集里的长度，随机取点就是在这个范围里面取
return len(self.label)
def __getitem__(self, i):
data=self.data[i]
label=self.label[i]
data = torch.tensor(data)
label = torch.tensor(label)
return data, label如果此文对你有帮助的话，还请不要吝惜你的点赞~
477
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
