Commit cbcac773 authored by Ngocson's avatar Ngocson

Fpaa simulation and capture

parent 2ee618ef
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 10 14:46:13 2019
@author: Ngocson
"""
import numpy as np
import matplotlib.pyplot as plt
import csv
global V
V = []
with open('ActivationFunction.csv') as csv_file:
csv_reader = csv.reader(csv_file, delimiter=';')
for row in csv_reader:
L = []
for value in row:
L.append(float(value))
V.append(L)
with open('santa_fe.txt') as sftr_file:
sftr_reader = csv.reader(sftr_file, delimiter=';')
santa_fe = np.array(next(sftr_reader)).astype(np.float)
santa_fe = santa_fe/santa_fe.max()
def first_order_numerical_filter(Y,X,tau = 5e-3):
return (X-Y)/tau
def approx_lin(x,vr=0.5,a=1.0,s=1.0):
global V
fi = min(max(20.0*vr,0.0),99.0)
if not fi.is_integer():
i = int(fi)
f = V[i]
else:
f = V[int(fi)]
fx = min(max(400.0*x*s,0.0),1999.0)
if not fx.is_integer():
x = int(fx)
return a*((x+1-fx)*(f[x]-f[x+1])+f[x+1])
else:
return a*f[int(fx)]
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_lin),
te = 10e-3,
tc = 10e-4,
N = 50,
vr = 0.3,
sp = None,
p = 0.1,
s_min = 0,
s_max = 1,
i_min = 0,
i_max = 1,
b_min = 0,
b_max = 1,
input_scalling=1):
s.activation = activation
s.f = f
s.te = te
s.tc = tc
s.N = N
s.sp = sp
s.p = p
s.input_scalling = input_scalling
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
if sp:
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)])
eigen,_ = np.linalg.eig(s.W)
s.W = s.W*s.sp/max(eigen.real)
else:
flag = True
while flag:
s.C = 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)])
s.W = np.matrix(np.zeros((s.N,s.N)))
flag = False
for i in range(s.N):
Ctot = s.C[i].sum()
if Ctot == 0:
flag = True
break
for j in range(s.N):
s.W[i,j] = 1.*s.C[i,j]/Ctot
plt.figure('weights')
plt.imshow(s.W)
plt.colorbar()
plt.show()
s.Ci = np.matrix([np.random.rand()*(s.i_max-s.i_min) + s.i_min for i in range(s.N)])
s.Wi = np.matrix(np.zeros(s.N))
Ctot = s.Ci.sum()
for i in range(s.N):
s.Wi[0,i] = 1.*s.Ci[0,i]/Ctot
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))))
s.nmos_voltages_dt = np.matrix(np.zeros(((1,s.N))))
s.filter_outputs_dt = np.matrix(np.zeros(((1,s.N))))
s.preactivation_dt = s.B
s.u_dt = 0
s.Vr = np.ones(s.N)*vr
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.preactivation_dt+(s.filter_outputs-s.filter_outputs_dt)*s.W+s.Wi*(u-s.u_dt)*s.input_scalling
s.u_dt = u
s.nmos_voltages_dt = s.nmos_voltages
s.filter_outputs_dt = s.filter_outputs
s.preactivation_dt = preactivation
s.nmos_voltages = s.activation(preactivation,s.Vr)
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,verbose = False):
outps = []
rep = []
for u in U:
s.step(u)
outps.append(s.read())
rep.append(s.filter_outputs[0,0])
if verbose:
plt.figure()
plt.title("Output")
plt.plot(rep)
plt.show()
return np.array(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)
s.Wo = tmp[:,:-1]
s.Bo = tmp[-1,-1]
X_ = (tmp*s.M)
_X = np.matrix(Xt)
plt.figure()
plt.title("State of the reservoir during training")
plt.imshow(s.M)
plt.colorbar()
plt.show()
plt.figure()
plt.title("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())
def response_NODE(s,U,n):
M = np.zeros(len(U))
for x in range(len(U)):
s.step(U[x])
M[x] =s.filter_outputs[0,n]
plt.figure()
plt.title("State of the node n°"+str(n))
plt.plot(M,'r-')
plt.show()
return M
def responses_NODE(s,U,N):
M = np.zeros((len(U),len(N)))
for x in range(len(U)):
s.step(U[x])
for n in range(len(N)):
M[x,n] =s.filter_outputs[0,N[n]]
plt.figure()
plt.title("State of the nodes "+str(N))
for n in range(len(N)):
plt.plot(M[:,n],label=str(N[n]))
plt.legend()
plt.show()
return M
N = 100
n_init = 2*N
n_tot = (2+7)*N
X = santa_fe[:-1]
Y = santa_fe[1:]
s = reservoirNMOS(N=N,
p = 5/N,
input_scalling=10.0)
###################### LET THE FUN BEGIN ##################
s.steps(X[:n_init])
e = s.fit(Y[n_init:n_tot],X[n_init:n_tot],0.001)
print(e)
s.reset(True)
Y_ = s.steps(X)
plt.figure()
plt.title('After training')
plt.plot(Y[n_init:],'r-',label='teacher')
plt.plot(Y_[n_init:],'b-',label='outp')
plt.legend()
plt.show()
plt.figure()
plt.title('Difference')
plt.plot(Y[n_init:]-Y[n_init-1:-1],'r-',label='difference on the laser')
plt.plot(Y[n_init:]-Y_[n_init:],'g-',label='difference on predicting')
plt.legend()
plt.show()
print(np.sqrt(np.linalg.norm(Y_[n_init:]-Y[n_init:])/Y_[n_init:].std()))
s.reset(True)
plt.figure()
plt.title('Testing')
plt.plot(Y[-100:],'r.-',label='teacher')
plt.plot(Y_[-100:],'b.-',label='outp')
plt.legend()
plt.show()
s.responses_NODE(X[:100],[1,2,3,4,5])
# 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
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
86;141;95;41;22;21;32;72;138;111;48;23;19;27;59;129;129;58;27;19;24;46;112;144;73;30;20;19;37;92;152;93;36;20;18;29;71;146;117;46;23;18;22;52;128;142;62;26;17;19;37;100;158;86;32;17;17;27;72;154;118;43;20;15;21;47;128;150;63;24;16;17;33;92;166;95;33;16;14;22;58;149;137;50;20;15;16;35;106;169;82;28;15;14;23;65;160;133;45;18;14;16;36;111;176;80;26;14;13;20;60;163;140;45;18;12;14;30;98;185;93;28;14;12;17;46;145;167;58;19;11;12;21;68;183;132;38;15;11;13;26;93;202;100;27;12;10;13;30;114;206;82;21;11;9;10;30;120;215;77;20;10;8;9;21;95;234;99;22;9;8;7;9;35;176;215;46;11;7;5;4;3;4;19;105;125;55;30;28;41;71;98;79;49;34;33;47;75;95;76;48;34;35;49;78;96;75;47;34;34;49;81;99;76;46;32;33;49;83;101;74;44;32;33;50;85;103;73;43;30;31;49;85;106;74;43;31;31;49;87;109;76;41;29;30;47;87;110;76;41;28;30;47;88;114;79;41;28;28;45;86;117;82;42;26;26;43;84;119;85;42;26;25;39;80;121;90;43;26;24;36;76;123;96;46;26;23;35;72;123;102;50;27;22;31;64;121;110;53;27;21;28;57;117;120;59;28;20;25;49;109;128;68;31;20;23;43;98;135;77;34;22;22;36;86;139;91;38;21;19;30;73;137;107;44;23;19;26;59;129;125;55;24;18;22;47;114;140;68;28;18;19;36;93;150;88;34;19;18;28;72;147;112;43;20;16;23;53;131;137;58;24;16;19;38;104;156;80;30;17;16;28;75;156;112;40;18;15;21;51;132;145;59;22;15;17;32;96;164;89;30;16;14;24;62;153;130;46;19;13;17;37;113;165;75;25;14;13;24;68;164;123;40;17;13;16;38;118;170;71;24;13;12;22;65;167;129;40;16;12;14;33;107;182;81;24;12;11;17;51;156;155;49;17;11;11;23;77;189;115;32;13;9;12;30;107;197;82;23;10;9;12;36;134;192;64;18;10;8;13;37;148;195;58;16;9;8;10;29;132;220;67;15;8;6;7;14;66;237;137;25;9;6;5;4;5;25;165;233;51;12;7;7;7;18;80;201;109;29;12;10;14;38;128;169;63;20;11;11;23;74;173;114;34;14;11;14;37;122;172;67;21;11;11;20;64;169;131;39;13;9;12;28;100;188;88;25;11;10;14;41;141;176;57;17;10;10;16;56;175;149;41;14;8;9;18;67;196;130;32;11;8;7;17;66;207;131;30;10;7;6;12;48;193;172;39;11;7;6;6;16;86;255;107;19;8;5;4;3;2;2;7;23;8;6;8;16;33;57;65;56;46;46;54;64;66;59;49;46;48;58;67;67;60;50;47;51;60;69;68;58;50;46;52;61;71;68;58;48;46;53;64;72;68;56;47;45;53;65;73;67;54;46;46;55;67;74;66;54;44;46;55;70;76;66;51;45;46;58;73;77;65;52;44;46;59;75;78;64;49;43;47;61;78;79;63;47;43;47;63;80;79;61;46;42;48;66;84;80;60;46;42;50;68;86;80;59;43;41;50;71;88;79;56;41;41;51;73;91;78;55;41;40;52;75;94;78;54;39;40;52;79;95;76;50;37;38;52;80;97;76;50;36;37;54;83;98;74;47;34;37;53;85;100;74;45;34;35;54;87;103;72;43;33;35;53;88;103;73;43;32;33;53;91;105;72;41;31;33;54;91;108;72;41;30;32;51;93;111;74;41;29;31;51;93;113;75;40;29;30;49;92;118;78;41;28;29;46;90;120;82;42;29;28;45;88;124;86;43;27;27;41;83;126;91;45;27;26;38;78;127;99;47;27;24;35;72;126;106;52;28;23;32;65;123;116;56;29;22;28;56;117;126;64;30;22;27;49;107;133;74;32;21;23;41;95;140;84;36;21;21;35;80;141;99;42;22;20;29;67;135;117;49;24;18;25;52;122;134;62;26;18;20;41;104;148;79;32;18;19;32;82;151;100;39;19;17;25;62;140;128;52;22;17;21;44;116;150;71;27;18;18;32;88;159;98;36;19;16;25;61;146;132;52;21;17;19;41;113;161;76;28;16;16;28;77;163;114;39;19;15;20;47;132;156;62;24;15;16;30;86;173;102;34;16;13;19;50;141;154;57;21;13;14;27;87;179;103;33;15;12;18;45;136;166;61;20;12;13;23
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Mon Jan 28 14:00:12 2019
@author: ngocson
"""
import numpy as np
class NARMA:
def __init__(s,n=10,alpha=0.3,beta=0.05,gamma=1.5,delta=0.1):
s.n = n
s.alpha = alpha
s.beta = beta
s.gamma = gamma
s.delta = delta
s.reset(np.zeros(n),(np.random.rand(n)-0.5)*0.2)
def reset(s,state,inputs):
s.state = state
s.inputs= inputs
def step(s,u):
yt = s.alpha*s.state[-1]\
+s.beta*s.state[-1]*s.state.sum()\
+s.gamma*s.inputs[0]*s.inputs[-1]\
+s.delta
s.inputs= np.append(s.inputs[1:],u)
s.state = np.append(s.state[1:],yt)
def steps(s,U):
for u in U:
s.step(u)
def read(s):
return s.state[-1]
def getTimeSerie(s,U):
Y = []
for u in U:
s.step(u)
Y.append(s.read())
return np.array(Y)
\ No newline at end of file
......@@ -6,12 +6,14 @@
import numpy as np
import matplotlib.pyplot as plt
import time
import csv
class RCnMOS:
def __init__(s,
dt = 10e-5,
wc = 5,
epsilon = 1e-4,
tol = 1e-9,
Vt = 0.5,
Kappa = 0.7,
Sigma = 0.05,
......@@ -38,7 +40,6 @@ class RCnMOS:
s.f = lambda U,t: (s.v0-U)*s.wc
s.tol = tol
s.epsilon = epsilon
......@@ -48,18 +49,10 @@ class RCnMOS:
def Vo(s,preactivation):
V = np.arange(s.vg,s.vdd,s.epsilon)
e = s.vdd - s.vg
mVo = 0
for vo in V:
I1 = s.I_vg_vd_vs(preactivation,s.vdd,vo)
I2 = s.I_vg_vd_vs(s.vr,vo,s.vg)
if abs(I1-I2) < s.tol:
return vo
elif abs(I1-I2) < e:
mVo = vo
e = abs(I1-I2)
return mVo
I1 = s.I_vg_vd_vs(preactivation,s.vdd,V)
I2 = s.I_vg_vd_vs(s.vr,V,s.vg)
dI = np.abs(I1-I2)
return V[np.argmin(dI)]
def activate(s,preactivation):
s.v0 = s.Vo(preactivation)
......@@ -79,25 +72,37 @@ class RCnMOS:
return s.u0
# FilterOutp, Activation, Previous Activation, preactvation
xo = np.array([0,0,0,0])
print("start")
Node = RCnMOS()
print("Node")
xX = []
plt.figure()
k=1
VR = [0,0.5,0.55,0.6,1,1.25,1.5,2,2.5]
VR = np.linspace(0,5,500)
for vr in VR:
X = []
Node.vr = vr
for i in range(500):
Node.activate(i/500*2.5)
print(k)
j = 0
start = time.time()
for u in np.linspace(0,5,5000):
j+=1
if j%500 == 0:
print(" u = ",u)
Node.activate(u)
X.append(Node.v0)
end = time.time()
xX.append(X)
plt.subplot(3,3,k)
plt.plot(np.linspace(0,2.5,500),X)
print("Calculation time: ",end-start)
plt.plot(np.linspace(0,5,5000),X,label="Vr = "+str(vr))
k+=1
plt.show()
with open('ActivationFunction.csv', mode='w') as file:
writer = csv.writer(file, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)
for L in xX:
writer.writerow(L)
'''
def van_der_pol(X,epsi,wo):
......
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 16 20:49:48 2019
@author: ngocson
"""
import numpy as np
from scipy import linalg
from random import random
import matplotlib.pyplot as plt
import csv
def write_matrix(matrix, filename = "matrice.csv"):
with open(filename, 'w') as csvfile:
writer = csv.writer(csvfile, delimiter=',')
for row in matrix:
writer.writerow(row)
class Reservoir:
def __init__(s,
N=100,
p=0.2,
sp=0.8,
outputScaling = 1,
v=0,
Activation = np.tanh,
Feedback = True,
inputScaling = 1,
Verbose = False,
damping = True):
s.N = N
s.v = v
s.p = p
s.sp = sp
s.outputScaling = outputScaling
s.inputScaling = inputScaling
s.Feedback = Feedback
s.Activation = Activation
s.initWeight()
if damping:
while not s.damped():
s.initWeight()
if Verbose:
print("W.std = ",s.W.std())
print("Wi.std = ",s.Wi.std())
ei,_ = np.linalg.eig(s.W)
print("W.eig.mean = ",np.abs(ei).mean())
print("W.eig.std = ",np.abs(ei).std())
ei,_ = np.linalg.eig(s.Wi.T*s.Wi)
print("Wi.eig.mean = ",np.abs(ei).mean())
print("Wi.eig.std = ",np.abs(ei).std())
if Feedback:
ei,_ = np.linalg.eig(s.Wfb.T*s.Wfb)
print("Wfb.eig.mean = ",np.abs(ei).mean())
print("Wfb.eig.std = ",np.abs(ei).std())
def initWeight(s):
#s.W = np.matrix([np.array( [0 if random()> p else -0.4 if random() < 0.5 else 0.4 for i in range(N)]) for j in range(s.N)])
s.W = np.matrix([np.array([random()-0.5 if random()<s.p else 0 for i in range(s.N)]) for j in range(s.N)])
#s.W = s.W - np.mean(s.W)
eigen,_ = linalg.eig(s.W)
s.W = s.W*s.sp/max(eigen.real)
s.reset()
s.Wi = np.matrix( np.random.rand(s.N) * 2 - 1)*s.inputScaling
#s.Wi = np.matrix([1 if random() < 0.5 else -1 for i in range(N)])
if s.Feedback:
s.Wfb = s.outputScaling*np.matrix(np.random.rand(s.N) * 2 - 1)
s.Wo = np.matrix([1.0 for i in range(s.N)])
s.Wio = np.matrix([1.0])
s.Bi = 1.0
def read(s,U):
return (s.Wo*s.state.T+s.Wio*U)[0,0] + s.Bi
def reset(s,Zeros = False):
#s.state = np.matrix(np.zeros((1,s.N)))
if Zeros:
s.state = np.matrix([0 for i in range(s.N)])
else:
s.state = np.matrix([2*random()-1 for i in range(s.N)])
def step(s,U=0,X = None): #U is a scalar
preactivation = s.state*s.W+ s.Wi*U
if s.Feedback:
if np.array_equal(X,np.array(None)):
preactivation += s.Wfb*s.read(U)
else:
preactivation += s.Wfb*X
s.state = s.Activation(preactivation) + s.v*(np.random.rand(s.N)-0.5)
return s.state
def steps(s,U, X = np.array(None)):
states = []
if np.array_equal(X,np.array(None)):
for u in U:
s.step(u)
states.append(s.read(u))
else:
for u,x in zip(U,X):
s.step(u,x)
states.append(s.read(u))
return states
def fit(s,Y,U,Forced = False, Show=False):
s.states_over_time = np.zeros((s.N+2,len(U)))
if s.Feedback:
s.step(U[0],Y[0])
s.states_over_time[:s.N+1,0] = np.concatenate((s.state.T,[[1.0]])).T
s.states_over_time[s.N+1:,0] = U[0]
for x in range(1,len(U)):
if Forced:
s.step(U[x],Y[x])
else:
s.step(U[x])
s.states_over_time[:s.N+1,x] = np.concatenate((s.state.T,[[1.0]])).T
s.states_over_time[s.N+1:,x] = U[x]
s.states_over_time = np.matrix(s.states_over_time.astype(float))
s.M = np.matrix(Y[1:])*np.linalg.pinv(s.states_over_time)
_X = np.matrix(Y[1:])
else:
for u in range(len(U)):
s.step(U[u])
s.states_over_time[:s.N+1,u] = np.concatenate((s.state.T,[[1.0]])).T
s.states_over_time[s.N+1:,u] = U[u]
s.states_over_time=np.matrix(s.states_over_time.astype(float))
s.M = np.matrix(Y)*np.linalg.pinv(s.states_over_time)
_X = np.matrix(Y)
s.Wo = s.M[0,:s.N]
s.Bi = s.M[0,s.N]
s.Wio = s.M[0,s.N+1:]
X_ = (s.M*s.states_over_time)
if Show:
plt.plot(_X.T)
plt.plot(X_.T)
plt.show()
return np.sqrt((np.array(X_-_X)**2).mean()/_X.std())
def damped(s):
s.reset()
M = np.zeros((s.N,s.N*2))
for i in range(s.N*2):
M[:,i] = s.state
s.step(0,0)
if ((M[:,-1]<0.01)*1).sum() < s.N:
return False
else:
if s.N <=50:
for k in range(1,1+s.N):
plt.subplot(int(np.sqrt(s.N)+1),int(np.sqrt(s.N)+1),k)
plt.plot(M[k-1,:])
plt.show()
return True
def start(s,U):
s.reset()
plt.plot(s.steps(U))
def evolution(s,U,X=np.array(None)):
s.reset()
M = np.zeros((s.N,len(U)))
for i in range(len(U)):
M[:,i] = s.state
if X != np.array(None):
s.step(U[i])
else:
s.step(U[i],X[i])
for k in range(1,1+s.N):
plt.subplot(int(np.sqrt(s.N)+1),int(np.sqrt(s.N)+1),k)
plt.plot(M[k-1,:])
plt.show()
def writeRes(s,fileprefix = "Reservoir_"):
write_matrix(s.W,fileprefix+"poidsInternes.csv")
write_matrix(s.Wo,fileprefix+"poidsSortie.csv")
write_matrix(s.Wi,fileprefix+"poidsEntree.csv")
write_matrix(s.Wfb,fileprefix+"poidsFeedback.csv")
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 16 15:20:31 2019