山东大学类脑实验 五 HMAX 模型实现
2024-04-21 10:54:34
山东大学类脑实验 五 HMAX 模型实现
实验目的：
加深对 HMAX 模型的理解，能够使用 HMAX 模型解决简单问题
实验内容：
根据 HMAX 模型的相关知识，使用 Python 语言实现一个简单的HMAX模型。
实验要求：
（1） 下载 MNIST 数据集。
（2） 构建 HMAX 模型。
（3） 使用 MNIST 数据集中的训练集训练网络，使用测试集测试训
练好的网络。
实验步骤：
所有代码文件在同一个目录下。
0.requirements
numpy
tqdm
torch
torchvision
sklearn
1.准备HMAX滤波器
下载HMAX滤波器参数文件（固定的，不需要训练），这里上传了一个免费的资源可以直接下载。https://download.csdn.net/download/ZEBRONE/85707155
然后创建hmax.py，内容为
# encoding: utf8
"""
PyTorch implementation of the HMAX model of human vision. For more information
about HMAX, check:http://maxlab.neuro.georgetown.edu/hmax.htmlThe S and C units of the HMAX model can almost be mapped directly onto
TorchVision's Conv2d and MaxPool2d layers, where channels are used to store the
filters for different orientations. However, HMAX also implements multiple
scales, which doesn't map nicely onto the existing TorchVision functionality.
Therefore, each scale has its own Conv2d layer, which are executed in parallel.Here is a schematic overview of the network architecture:layers consisting of units with increasing scale
S1 S1 S1 S1 S1 S1 S1 S1 S1 S1 S1 S1 S1 S1 S1 S1\ / \ / \ / \ / \ / \ / \ / \ /C1 C1 C1 C1 C1 C1 C1 C1\ \ \ | / / / /ALL-TO-ALL CONNECTIVITY/ / / | \ \ \ \S2 S2 S2 S2 S2 S2 S2 S2| | | | | | | |C2 C2 C2 C2 C2 C2 C2 C2Author: Marijn van Vliet <w.m.vanvliet@gmail.com>References
----------.. [1] Riesenhuber, Maximilian, and Tomaso Poggio. “Hierarchical Models ofObject Recognition in Cortex.” Nature Neuroscience 2, no. 11 (1999):1019–25. https://doi.org/10.1038/14819... [2] Serre, T, M Kouh, C Cadieu, U Knoblich, Gabriel Kreiman, and T Poggio.“A Theory of Object Recognition: Computations and Circuits in theFeedforward Path of the Ventral Stream in Primate Visual Cortex.”Artificial Intelligence, no. December (2005): 1–130.https://doi.org/10.1.1.207.9279... [3] Serre, Thomas, Aude Oliva, and Tomaso Poggio. “A FeedforwardArchitecture Accounts for Rapid Categorization.” Proceedings of theNational Academy of Sciences 104, no. 15 (April 10, 2007): 6424–29.https://doi.org/10.1073/pnas.0700622104... [4] Serre, Thomas, and Maximilian Riesenhuber. “Realistic Modeling ofSimple and Complex Cell Tuning in the HMAXModel, and Implications forInvariant Object Recognition in Cortex.” CBCL Memo, no. 239 (2004)... [5] Serre, Thomas, Lior Wolf, Stanley Bileschi, Maximilian Riesenhuber,and Tomaso Poggio. “Robust Object Recognition with Cortex-likeMechanisms.” IEEE Trans Pattern Anal Mach Intell 29, no. 3 (2007):411–26. https://doi.org/10.1109/TPAMI.2007.56.
"""
import numpy as np
from scipy.io import loadmat
import torch
from torch import nndef gabor_filter(size, wavelength, orientation):"""Create a single gabor filter.Parameters----------size : intThe size of the filter, measured in pixels. The filter is square, henceonly a single number (either width or height) needs to be specified.wavelength : floatThe wavelength of the grating in the filter, relative to the half thesize of the filter. For example, a wavelength of 2 will generate aGabor filter with a grating that contains exactly one wave. Thisdetermines the "tightness" of the filter.orientation : floatThe orientation of the grating in the filter, in degrees.Returns-------filt : ndarray, shape (size, size)The filter weights."""lambda_ = size * 2. / wavelengthsigma = lambda_ * 0.8gamma = 0.3 # spatial aspect ratio: 0.23 < gamma < 0.92theta = np.deg2rad(orientation + 90)# Generate Gabor filterx, y = np.mgrid[:size, :size] - (size // 2)rotx = x * np.cos(theta) + y * np.sin(theta)roty = -x * np.sin(theta) + y * np.cos(theta)filt = np.exp(-(rotx**2 + gamma**2 * roty**2) / (2 * sigma ** 2))filt *= np.cos(2 * np.pi * rotx / lambda_)filt[np.sqrt(x**2 + y**2) > (size / 2)] = 0# Normalize the filterfilt = filt - np.mean(filt)filt = filt / np.sqrt(np.sum(filt ** 2))return filtclass S1(nn.Module):"""A layer of S1 units with different orientations but the same scale.The S1 units are at the bottom of the network. They are exposed to the rawpixel data of the image. Each S1 unit is a Gabor filter, which detectsedges in a certain orientation. They are implemented as PyTorch Conv2dmodules, where each channel is loaded with a Gabor filter in a specificorientation.Parameters----------size : intThe size of the filters, measured in pixels. The filters are square,hence only a single number (either width or height) needs to bespecified.wavelength : floatThe wavelength of the grating in the filter, relative to the half thesize of the filter. For example, a wavelength of 2 will generate aGabor filter with a grating that contains exactly one wave. Thisdetermines the "tightness" of the filter.orientations : list of floatThe orientations of the Gabor filters, in degrees."""def __init__(self, size, wavelength, orientations=[90, -45, 0, 45]):super().__init__()self.num_orientations = len(orientations)self.size = size# Use PyTorch's Conv2d as a base object. Each "channel" will be an# orientation.self.gabor = nn.Conv2d(1, self.num_orientations, size,padding=size // 2, bias=False)# Fill the Conv2d filter weights with Gabor kernels: one for each# orientationfor channel, orientation in enumerate(orientations):self.gabor.weight.data[channel, 0] = torch.Tensor(gabor_filter(size, wavelength, orientation))# A convolution layer filled with ones. This is used to normalize the# result in the forward method.self.uniform = nn.Conv2d(1, 4, size, padding=size // 2, bias=False)nn.init.constant_(self.uniform.weight, 1)# Since everything is pre-computed, no gradient is requiredfor p in self.parameters():p.requires_grad = Falsedef forward(self, img):"""Apply Gabor filters, take absolute value, and normalize."""s1_output = torch.abs(self.gabor(img))norm = torch.sqrt(self.uniform(img ** 2))norm.data[norm == 0] = 1 # To avoid divide by zeros1_output /= normreturn s1_outputclass C1(nn.Module):"""A layer of C1 units with different orientations but the same scale.Each C1 unit pools over the S1 units that are assigned to it.Parameters----------size : intSize of the MaxPool2d operation being performed by this C1 layer."""def __init__(self, size):super().__init__()self.size = sizeself.local_pool = nn.MaxPool2d(size, stride=size // 2,padding=size // 2)def forward(self, s1_outputs):"""Max over scales, followed by a MaxPool2d operation."""s1_outputs = torch.cat([out.unsqueeze(0) for out in s1_outputs], 0)# Pool over all scaless1_output, _ = torch.max(s1_outputs, dim=0)# Pool over local (c1_space x c1_space) neighbourhoodreturn self.local_pool(s1_output)class S2(nn.Module):"""A layer of S2 units with different orientations but the same scale.The activation of these units is computed by taking the distance betweenthe output of the C layer below and a set of predefined patches. Thisdistance is computed as:d = sqrt( (w - p)^2 )= sqrt( w^2 - 2pw + p^2 )Parameters----------patches : ndarray, shape (n_patches, n_orientations, size, size)The precomputed patches to lead into the weights of this layer.activation : 'gaussian' | 'euclidean'Which activation function to use for the units. In the PNAS paper, agaussian curve is used ('guassian', the default), whereas the MATLABimplementation of The Laboratory for Computational CognitiveNeuroscience uses the euclidean distance ('euclidean').sigma : floatThe sharpness of the tuning (sigma in eqn 1 of [1]_). Defaults to 1.References:-----------.. [1] Serre, Thomas, Aude Oliva, and Tomaso Poggio. “A FeedforwardArchitecture Accounts for Rapid Categorization.” Proceedings of theNational Academy of Sciences 104, no. 15 (April 10, 2007): 6424–29.https://doi.org/10.1073/pnas.0700622104."""def __init__(self, patches, activation='gaussian', sigma=1):super().__init__()self.activation = activationself.sigma = sigmanum_patches, num_orientations, size, _ = patches.shape# Main convolution layerself.conv = nn.Conv2d(in_channels=num_orientations,out_channels=num_orientations * num_patches,kernel_size=size,padding=size // 2,groups=num_orientations,bias=False)self.conv.weight.data = torch.Tensor(patches.transpose(1, 0, 2, 3).reshape(1600, 1, size, size))# A convolution layer filled with ones. This is used for the distance# computationself.uniform = nn.Conv2d(1, 1, size, padding=size // 2, bias=False)nn.init.constant_(self.uniform.weight, 1)# This is also used for the distance computationself.patches_sum_sq = nn.Parameter(torch.Tensor((patches ** 2).sum(axis=(1, 2, 3))))self.num_patches = num_patchesself.num_orientations = num_orientationsself.size = size# No gradient required for this layerfor p in self.parameters():p.requires_grad = Falsedef forward(self, c1_outputs):s2_outputs = []for c1_output in c1_outputs:conv_output = self.conv(c1_output)# Unstack the orientationsconv_output_size = conv_output.shape[3]conv_output = conv_output.view(-1, self.num_orientations, self.num_patches, conv_output_size,conv_output_size)# Pool over orientationsconv_output = conv_output.sum(dim=1)# Compute distancec1_sq = self.uniform(torch.sum(c1_output ** 2, dim=1, keepdim=True))dist = c1_sq - 2 * conv_outputdist += self.patches_sum_sq[None, :, None, None]# Apply activation functionif self.activation == 'gaussian':dist = torch.exp(- 1 / (2 * self.sigma ** 2) * dist)elif self.activation == 'euclidean':dist[dist < 0] = 0 # Negative values should never occurtorch.sqrt_(dist)dist = -distelse:raise ValueError("activation parameter should be either ""'gaussian' or 'euclidean'.")s2_outputs.append(dist)return s2_outputsclass C2(nn.Module):"""A layer of C2 units operating on a layer of S2 units."""def forward(self, s2_outputs):"""Take the maximum value of the underlying S2 units."""maxs = [s2.max(dim=3)[0] for s2 in s2_outputs]maxs = [m.max(dim=2)[0] for m in maxs]maxs = torch.cat([m[:, None, :] for m in maxs], 1)# maxs = torch.where(torch.isnan(maxs), torch.full_like(maxs, -100), maxs)return maxs.max(dim=1)[0]class HMAX(nn.Module):"""The full HMAX model.Use the `get_all_layers` method to obtain the activations for all layers.If you are only interested in the final output (=C2 layer), use the modelas any other PyTorch module:model = HMAX(universal_patch_set)output = model(img)Parameters----------universal_patch_set : strFilename of the .mat file containing the universal patch set.s2_act : 'gaussian' | 'euclidean'The activation function for the S2 units. Defaults to 'gaussian'.Returns-------c2_output : list of Tensors, shape (batch_size, num_patches)For each scale, the output of the C2 units."""def __init__(self, universal_patch_set, s2_act='gaussian'):super().__init__()# S1 layers, consisting of units with increasing sizeself.s1_units = [S1(size=7, wavelength=4),S1(size=9, wavelength=3.95),S1(size=11, wavelength=3.9),S1(size=13, wavelength=3.85),S1(size=15, wavelength=3.8),S1(size=17, wavelength=3.75),S1(size=19, wavelength=3.7),S1(size=21, wavelength=3.65),S1(size=23, wavelength=3.6),S1(size=25, wavelength=3.55),S1(size=27, wavelength=3.5),S1(size=29, wavelength=3.45),S1(size=31, wavelength=3.4),S1(size=33, wavelength=3.35),S1(size=35, wavelength=3.3),S1(size=37, wavelength=3.25),]# Explicitly add the S1 units as submodules of the modelfor s1 in self.s1_units:self.add_module('s1_%02d' % s1.size, s1)# Each C1 layer pools across two S1 layersself.c1_units = [C1(size=8),C1(size=10),C1(size=12),C1(size=14),C1(size=16),C1(size=18),C1(size=20),C1(size=22),]# Explicitly add the C1 units as submodules of the modelfor c1 in self.c1_units:self.add_module('c1_%02d' % c1.size, c1)# Read the universal patch set for the S2 layerm = loadmat(universal_patch_set)patches = [patch.reshape(shape[[2, 1, 0, 3]]).transpose(3, 0, 2, 1)for patch, shape in zip(m['patches'][0], m['patchSizes'].T)]# One S2 layer for each patch scale, operating on all C1 layersself.s2_units = [S2(patches=scale_patches, activation=s2_act)for scale_patches in patches]# Explicitly add the S2 units as submodules of the modelfor i, s2 in enumerate(self.s2_units):self.add_module('s2_%d' % i, s2)# One C2 layer operating on each scaleself.c2_units = [C2() for s2 in self.s2_units]# Explicitly add the C2 units as submodules of the modelfor i, c2 in enumerate(self.c2_units):self.add_module('c2_%d' % i, c2)def run_all_layers(self, img):"""Compute the activation for each layer.Parameters----------img : Tensor, shape (batch_size, 1, height, width)A batch of images to run through the modelReturns-------s1_outputs : List of Tensors, shape (batch_size, num_orientations, height, width)For each scale, the output of the layer of S1 units.c1_outputs : List of Tensors, shape (batch_size, num_orientations, height, width)For each scale, the output of the layer of C1 units.s2_outputs : List of lists of Tensors, shape (batch_size, num_patches, height, width)For each C1 scale and each patch scale, the output of the layer ofS2 units.c2_outputs : List of Tensors, shape (batch_size, num_patches)For each patch scale, the output of the layer of C2 units."""s1_outputs = [s1(img) for s1 in self.s1_units]# Each C1 layer pools across two S1 layersc1_outputs = []for c1, i in zip(self.c1_units, range(0, len(self.s1_units), 2)):c1_outputs.append(c1(s1_outputs[i:i+2]))s2_outputs = [s2(c1_outputs) for s2 in self.s2_units]c2_outputs = [c2(s2) for c2, s2 in zip(self.c2_units, s2_outputs)]return s1_outputs, c1_outputs, s2_outputs, c2_outputsdef forward(self, img):"""Run through everything and concatenate the output of the C2s."""c2_outputs = self.run_all_layers(img)[-1]c2_outputs = torch.cat([c2_out[:, None, :] for c2_out in c2_outputs], 1)return c2_outputsdef get_all_layers(self, img):"""Get the activation for all layers as NumPy arrays.Parameters----------img : Tensor, shape (batch_size, 1, height, width)A batch of images to run through the modelReturns-------s1_outputs : List of arrays, shape (batch_size, num_orientations, height, width)For each scale, the output of the layer of S1 units.c1_outputs : List of arrays, shape (batch_size, num_orientations, height, width)For each scale, the output of the layer of C1 units.s2_outputs : List of lists of arrays, shape (batch_size, num_patches, height, width)For each C1 scale and each patch scale, the output of the layer ofS2 units.c2_outputs : List of arrays, shape (batch_size, num_patches)For each patch scale, the output of the layer of C2 units."""s1_out, c1_out, s2_out, c2_out = self.run_all_layers(img)return ([s1.cpu().detach().numpy() for s1 in s1_out],[c1.cpu().detach().numpy() for c1 in c1_out],[[s2_.cpu().detach().numpy() for s2_ in s2] for s2 in s2_out],[c2.cpu().detach().numpy() for c2 in c2_out],)2.下载MNIST数据并用HMAX滤波器处理数据(GTX1650 用时35min）
import torch
import torchvision
from torch.utils.data import DataLoader
from torchvision import datasets, transforms
import pickle
import numpy as np
import hmax
from tqdm import tqdm# Initialize the model with the universal patch set
print('Constructing model')
model = hmax.HMAX('./universal_patch_set.mat')path = './data/'
transform = torchvision.transforms.Compose([torchvision.transforms.Grayscale(),torchvision.transforms.ToTensor(),torchvision.transforms.Lambda(lambda x: x * 255)])
# torchvision.transforms.Normalize(mean=[0.5], std=[0.5])])trainData = torchvision.datasets.MNIST(path, train=True, transform=transform, download=True)
testData = torchvision.datasets.MNIST(path, train=False, transform=transform)BATCH_SIZE = 1
trainDataLoader = torch.utils.data.DataLoader(dataset=trainData, batch_size=BATCH_SIZE, shuffle=True)
testDataLoader = torch.utils.data.DataLoader(dataset=testData, batch_size=BATCH_SIZE)# Determine whether there is a compatible GPU available
device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')print('Running model on', device)model = model.to(device)trainData_after = []
for X, y in tqdm(trainDataLoader):c2 = model.forward(X.to(device))[:, 0, :]trainData_after.append((c2[0].cpu(), y))
np.save(open("./data/train_data_new.npy", "wb"), trainData_after)
testData_after = []
for X, y in tqdm(testDataLoader):c2 = model.forward(X.to(device))[:, 0, :]testData_after.append((c2[0].cpu(), y))
np.save(open("./data/test_data_new.npy", "wb"), testData_after)
3.用SVM作为VTU训练（GTX1650 用时10min）
from sklearn.svm import SVC
import numpy as npprint("loading data...")
trainData = np.load(open("./data/train_data_new.npy","rb"),allow_pickle=True)
testData = np.load(open("./data/test_data_new.npy","rb"),allow_pickle=True)
print("converting data...")
X = [_[0].cpu().numpy().tolist() for _ in trainData]
y = [int(_[1]) for _ in trainData]
X_test = [_[0].cpu().numpy().tolist() for _ in testData]
y_test = [int(_[1]) for _ in testData]print("training data...")
model = SVC(kernel="rbf")
model.fit(X,y)print("testing acc:",model.score(X_test,y_test))4.训练完成直接显示测试结果
testing acc: 0.9132（恒定值，因为是SVM）
山东大学类脑实验 五 HMAX 模型实现相关推荐
- 【类脑实验】`Hopfield` 模型的实现
类脑实验记录系列:实验1 Hopfield 模型的实现 实验名称:Hopfield 模型的实现 课程名称:认知科学与类脑计算 一 实验目的: 加深对 Hopfield 模型的理解,能够使用 Hopfi ...
- 山东大学linux应用实验五,【Linux】山东大学Linux应用课程实验记录
找到这篇博文的人,一定被Linux实验弄得很爆炸吧哈哈哈. 这里是我Linux实验的记录,供大家学习和参考.如有错误,还请指正. 实验一 一. 基本命令 显示系统当前时间. date 显示2003年的 ...
- matlab复杂周期信号类建立,实验五 基于Matlab的信号频谱分析（复杂）
本次实验注意:<实验五MALTAB基础知识(简单)> <实验五 基于Matlab的信号频谱分析(复杂)> 选作一个即可 实验五 基于Matlab的信号频谱分析 (一) 实验目的 ...
- 山东大学软件学院数据挖掘实验五（2）的坑
一. 实验目的 掌握数据导入Hive表的方式 理解三种数据导入Hive表的原理 二. 实验内容 1.启动Hadoop和Hive服务并创建数据表 2.将Hive表中的数据导出 三. 实验 ...
- 山东大学软件学院数据库系统实验五
文章目录 一.实验时间 二.实验题目 一.实验时间 2021年5月4日星期二,第10周 二.实验题目 1． 在学生表pub.student中统计名字(姓名的第一位是姓氏,其余为名字,不考虑复姓)的使用 ...
- 山东大学Linux应用实验五
本文属于原创,转载请注明出处. 实验步骤与内容: 一.引号的使用 双引号和单引号的使用.依次输入下列命令: (1) s t r i n g = " h e l l o w o r l d ! ...
- 清华“天机芯”团队再发重磅研究！以全新类脑计算系统实现通用人工智能
关注ITValue,看企业级最新鲜.最价值报道! 10月15日,清华大学计算机系张悠慧团队和精仪系施路平团队与合作者发表一项最新类脑计算体系结构的突破性研究成果,首次提出"类脑计算完备性&q ...
- 清华团队类脑计算首推「神经形态完备性」，通用人工智能要来了
2020-10-16 14:13:46 [新智元导读]人脑,是自然界中最完美的计算系统,不仅能处理多种复杂任务还有最高的效能,是目前唯一的「通用智能体」.最近,清华联合团队提出的类脑计算系统新框架,或 ...
- 清华类脑计算成果再登Nature：张悠慧施路平团队出品，有望打破冯诺依曼瓶颈...
贾浩楠 发自 凹非寺 量子位 报道 | 公众号 QbitAI 清华类脑计算研究成果,再登Nature. 新研究的关键词是:类脑计算.新计算机系统框架.通用人工智能(AGI). 它的重要性,在于有希望打 ...
最新文章
- H3 BPM钉钉接入配置
- R语言使用caret包构建gbdt模型（随机梯度提升树、Stochastic Gradient Boosting ）构建回归模型、通过method参数指定算法名称
- 人物访谈：松本行弘谈Ruby
- WeChat报错微信小程序图片加载失败渲染层网络层错误Failed to load image /pages/index/image/index.jpg：用绝对路径不用相对路径
- OpenSitUp开源项目：零基础开发基于姿态估计的运动健身APP
- 使用PowerShell配置Microsoft Teams
- Dnslog在SQL注入中的利用
- (Android小应用)在Android中实现多线程断点下载（连载二）
- python替换文件中的字符串_Python文件操作中进行字符串替换的方法（保存到新文件/当前文件）...
- pyspark分类算法之决策树分类器模型实践【decisionTreeClassifier】
- 6-5.添加HLSL顶点着色
- 数据中心运维管理方案
- VS2019离线安装方法
- DONET中常用的一些快捷键收集。
- 游戏公司可能碰到哪些知识产权问题?
- linux mysql.sock文件_关于linux上mysql.sock文件的个人理解
- Choerodon猪齿鱼实践之持续交付流水线
- 【Java】命名规范
- Git常用命令及方法大全
- python学习 -对象把微信消息撤回后好慌，有了这个你就能看到撤回的消息了（超详解）