-
Notifications
You must be signed in to change notification settings - Fork 2
/
dnf_1d.py
202 lines (161 loc) · 6.66 KB
/
dnf_1d.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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
##########################################################
## ##
## Module: dnf_1d.py ##
## ##
## Version: 0.2 ##
## ##
## Description: A running simulation of a dynamic ##
## neural field with several modifications for ##
## improving runtime. This is the 1d version. It ##
## is used as a simple starting point for underst- ##
## ing path integration along with a working ##
## implementation of the mexican hat gaussians which ##
## enable multiple hypothesese on a single neural ##
## field ##
## ##
## Checklist: 1- Enable multi bubbles ##
## 2- Enable path integration ##
## ##
##########################################################
from numpy import *
from scipy.signal import *
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib.ticker import LinearLocator, FormatStrFormatter
class dnf:
######### Simulation variables ############
n=100 #Number of nodes per side
tau=0.1 #Tau
I=zeros((n,1)) #Input activity
dx=2*pi/n #Length per node in the x dimension
sig=2*pi/10*0.6 #Sigma of the gaussian function
c= 0.095 #Global inhibition
nrate = 0.3
###########################################
#Takes a location between 0 and 2pi for both x along with the sigma
# value for the gaussian function. Produces an narray with the gaussian
# figure.
@staticmethod
def gauss_pbc(locx,sig):
w=zeros((dnf.n,1))
for i in range(dnf.n):
d=min([abs(i*dnf.dx-locx) , 2*pi-abs(i*dnf.dx-locx)])
w[i]=1./(sqrt(2*pi)*sig)\
*exp(-(d**2/(2*sig**2)))
return w
@staticmethod
def gauss_diff(locx,sig1,sig2):
w=zeros((dnf.n,1))
for i in range(dnf.n):
d=min([abs(i*dnf.dx-locx) , 2*pi-abs(i*dnf.dx-locx)])
w[i]=1./(sqrt(2*pi)*sig1)*exp(-(d**2/(2*sig1**2)))\
- 1./(sqrt(2*pi)*sig2)*exp(-(d**2/(2*sig2**2)))
return w
#Initializor for the dynamic neural field. It creates a new one with the
# parameters described above. In later versions it would support a multi
# hypothesis dnf by supplying it with the parameter multi=True
def __init__(self,multi=False):
self.t = 0
self.gauss = dnf.gauss_pbc(pi,self.sig)
self.u=zeros((dnf.n,1)) #The neural field state at a specific time
if multi:
self.c = 0.03
self.z = 1000*(self.hebbMulti()-self.c)
else:
self.z = 1000*(self.hebb()-self.c)
self.zpi1 = 1000*(self.hebbPI1()-self.c)
self.zpi2 = 1000*(self.hebbPI2()-self.c)
self.xall=[]
#Produces the weight array that is responsible for integrating motion
# to the right
def hebbPI1(self):
w=zeros((self.n,self.n))
rbar = dnf.gauss_pbc(0,self.sig)
for i in range(self.n):
r=dnf.gauss_pbc(i*self.dx,self.sig)
w=w+dot(r,rbar.transpose())
rbar= (1-self.nrate)*rbar+ self.nrate*r
return w/self.n
#Produces the weight array that is responsible for integrating motion
# to the left
def hebbPI2(self):
w=zeros((self.n,self.n))
rbar = dnf.gauss_pbc(0,self.sig)
for i in range(self.n):
r=dnf.gauss_pbc((self.n-i)*self.dx,self.sig)
w=w+dot(r,rbar.transpose())
rbar= (1-self.nrate)*rbar+ self.nrate*r
return w/self.n
#Uses hebbian learning to produce a single n*n array of the relation between
# two gaussians. It yields the weights between the middle node and all the
# others.
def hebb(self):
w=zeros((self.n,self.n))
for i in range(self.n):
r=dnf.gauss_pbc(i*self.dx,self.sig)
w=w+dot(r,r.transpose())
return w/self.n
#Produces a weight array which supports multi hypotheses.
def hebbMulti(self):
w=zeros((self.n,self.n))
for i in range(self.n):
r=dnf.gauss_pbc(i*self.dx,self.sig/3) - dnf.gauss_pbc(i*self.dx,self.sig)
w=w+dot(r,r.transpose())
return w/self.n
#Takes a dynamic field state along with an activity input and the weight
# from hebb. It updates the dynamic neural field and returns it after one
# step in time.
def update(self,I, v=0):
self.t += 1
if v <= 0:
zpi = self.zpi2
else:
zpi = self.zpi1
r=0.5*(tanh(0.1*self.u)+1)
z = (self.z + abs(v)*zpi)/(abs(v)+1)
self.u=self.u+ self.tau*(-self.u+dot(z,r)*self.dx+I)
x=0.5*(tanh(0.1*self.u)+1)
self.xall = self.xall + [list(self.u.reshape(-1,))]
#Takes a dnf state and draws it in 3dwith respect to time.
def plot(self):
fig = plt.figure()
ax = fig.gca(projection='3d')
self.X, self.Y = meshgrid(
arange(self.n),
arange(self.t)) #A grid for plotting purposes
surf = ax.plot_surface(self.X, self.Y, self.xall, cmap=cm.jet,
linewidth=0, antialiased=True)
ax.w_zaxis.set_major_locator(LinearLocator(10))
ax.w_zaxis.set_major_formatter(FormatStrFormatter('%.1f'))
ax.set_zlabel("Excitation")
ax.set_xlabel("Node 'X'")
ax.set_ylabel("Time")
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
#Test case for the module
if __name__ == "__main__":
dnfex = dnf(multi=True)
# r = dnf.gauss_pbc(pi,dnf.sig/2) - dnf.gauss_pbc(pi,dnf.sig)
# plt.plot(r)
# plt.show()
#Provide input at pi/2 for 50 steps
I=dnf.gauss_pbc(pi,dnf.sig/10)
for i in range(200):
dnfex.update(I)
I=zeros((dnf.n,1))
for i in range(200):
dnfex.update(I)
I=dnf.gauss_pbc(2*pi/3,dnf.sig/10)
for i in range(300):
dnfex.update(I)
I=zeros((dnf.n,1))
for i in range(20):
dnfex.update(I)
I=dnf.gauss_pbc(2*pi,dnf.sig/20)
for i in range(200):
dnfex.update(I)
I=zeros((dnf.n,1))
for i in range(200):
dnfex.update(I)
dnfex.plot()