...
 
import numpy as np
class LeakyEsn:
def __init__(s,
N,
Sp,
p,
a,
c,
dt,
Activation,
v,
inputScalling,
outputScalling,
Wmin,
Wmax):
s.N = N
s.Activation = Activation
s.v = v
s.p = p
s.sp = sp
s.a = a
s.c = c
s.dt = dt
s.outputScalling= outputScalling
s.inputScalling = inputScalling
s.Wmin = Wmin
s.Wmax = Wmax
s.states = np.zeros(s.N)
s.W = Ïnp.matrix( [ [(s.Wmax-s.Wmin)*np.random.rand() + s.Wmin for i in range(s.N)] for j in range(s.N)] )
while not s.damped(1e-4):
s.W = np.matrix( [ [(s.Wmax-s.Wmin)*np.random.rand() + s.Wmin for i in range(s.N)] for j in range(s.N)] )
pass
def damped(s,e)
s.states = np.Random.rand(s.N)
Y = steps(np.zeros((s.N,10*s.N)))
if np.sum((Y[-1]>e).astype(np.int))>0:
return False
return True
def steps(s,U,X):
pass
def step(s,U):
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 10 14:46:13 2019
@author: Ngocson
"""
import numpy as np
import matplotlib.pyplot as plt
def first_order_numerical_filter(Y,X,tau = 1e-3):
return (X-Y)/tau
def Approx(x):
P = [-0.0028984 , 0.02966824, -0.11060687, 0.27304077, -0.1222173 , 0.01029023]
def PolyEval(x,P):
y = 0
for i in range(len(P)):
y += x**(len(P)-i-1)*P[i]
return y
if x < 0.11:
return 0
else:
return min(2.5,PolyEval(x,P))
class reservoirNMOS:
def __init__(s ,
f = first_order_numerical_filter,
activation = np.vectorize(Approx),
te = 10e-3,
tc = 10e-5,
N = 50,
sp = None,
p = 0.1,
s_min = 0,
s_max = 1,
i_min = 0,
i_max = 1,
b_min = 0,
b_max = 1,
n_states = 2):
s.activation = activation
s.f = f
s.te = te
s.tc = tc
s.N = N
s.n_states = n_states
s.sp = sp
s.p = p
s.s_min = s_min
s.s_max = s_max
s.i_max = i_max
s.i_min = i_min
s.b_min = b_min
s.b_max = b_max
s.W = np.matrix([[np.random.rand()*(s.s_max-s.s_min) + s.s_min if np.random.rand() < p else 0 for j in range(s.N)] for i in range(s.N)])
if sp:
eigen,_ = np.linalg.eig(s.W)
s.W = s.W*s.sp/max(eigen.real)
plt.figure('weights')
plt.imshow(s.W)
plt.colorbar()
plt.show()
s.Wi = np.matrix([np.random.rand()*(s.i_max-s.i_min) + s.i_min for i in range(s.N)])
s.B = np.matrix([np.random.rand()*(s.b_max-s.b_min) + s.b_min for i in range(s.N)])
s.Wo = np.matrix(np.zeros((1,s.N)))
s.Bo = 0
s.nmos_voltages = np.matrix(np.zeros(((1,s.N))))
s.filter_outputs = np.matrix(np.zeros(((1,s.N))))
def reset(s,rand = False):
if rand:
s.filter_outputs = np.matrix(np.random.random(((1,s.N))))
else:
s.filter_outputs = np.matrix(np.zeros(((1,s.N))))
def read(s):
return (s.Wo*s.filter_outputs.T)[0,0] + s.Bo
def m_step(s,u):
preactivation = s.filter_outputs*s.W+s.B+s.Wi*u
s.nmos_voltages = s.activation(preactivation)
k1 = s.f(s.filter_outputs,
s.nmos_voltages)
k2 = s.f(s.filter_outputs+0.5*k1*s.tc,
s.nmos_voltages)
k3 = s.f(s.filter_outputs+0.5*k2*s.tc,
s.nmos_voltages)
k4 = s.f(s.filter_outputs+k1*s.tc,
s.nmos_voltages)
s.filter_outputs = s.filter_outputs + s.tc/6.*(k1 + 2*k2 + 2*k3 + k4)
def step(s,u):
n_iter = s.te/s.tc
for i in range(int(n_iter)):
s.m_step(u)
def steps(s,U):
outps = []
rep = []
for u in U:
s.step(u)
outps.append(s.read())
rep.append(s.filter_outputs[0,0])
print(rep)
plt.figure()
plt.plot(rep)
plt.show()
return outps
def fit(s,Xt,U,lmbda = 0.05):
s.M = np.zeros((s.N+1,len(U)))
for x in range(len(U)):
s.step(U[x])
s.M[:,x] = np.concatenate((s.filter_outputs.T,[[1]])).T
s.M = np.matrix(s.M.astype(float))
tmp = np.matrix(Xt)*(np.linalg.inv(s.M.T*s.M+lmbda*np.eye(Xt.shape[0]))*s.M.T)
print(tmp.shape)
s.Wo = tmp[:,:-1]
s.Bo = tmp[-1,-1]
X_ = (tmp*s.M)
_X = np.matrix(Xt)
plt.figure("training results")
plt.plot(_X.T,'r-',label="teacher")
plt.plot(X_.T,'b-', label="reservoir output")
plt.legend()
plt.show()
return np.sqrt(np.linalg.norm(X_-_X)/X_.var())
N = 50
n_init = 2*N
n_tot = 10*N+n_init
X= (np.sin(np.linspace(1,12,n_tot))+1)*0.25
Y = (X-1)**7+1
s = reservoirNMOS(N=N)
s.steps(X[:n_init])
e = s.fit(Y[n_init:],X[n_init:],0.001)
s.reset(True)
# s.steps(X[:n_init])
# Fonction différentes non-linéarité: rentrer un tableau de valeur avec interpollation linéaire
# Améliorer notre compréhension du système:
# Afficher la réponse de noeuds excitation
# Doubler le temps d'entrée (bloqué sur deux périodes)
# Quantité de noeuds
# Input scalling
#Test
# Santafe database -> prédiction de série chaotiques à un pas : NMSE ~10^-2
# Fonction on linéaires aléatoires
X= (np.sin(np.linspace(1,36,3*n_tot))+1)
Y = (X-1)**7+1
Y_ = s.steps(X)
plt.figure()
plt.plot(Y_,'b.-',label='outp')
plt.plot(Y,'r.-',label='teacher')
plt.show()
print(e)
\ No newline at end of file