Commit f8820280 authored by Sayah-Sel's avatar Sayah-Sel
Browse files

several epsilon using autograd

parent 86359a34
import numpy as np
import autograd.numpy as npa
import scipy.spatial as sp
from scipy.linalg import cholesky, cho_solve
from numpy import fft as f
from autograd import grad, jacobian
from autograd import elementwise_grad as egrad
from tqdm import tqdm
import multiprocessing as mp
from tqdm.contrib.concurrent import process_map
import os
d=2
#dx=0.02
L=1
N=32
dx=L/N
#N=L/dx
mu=0
corr_L=["norm1","norm2"]
corr="norm1"
l_cV=0.25
l_cA=0.25
#e_L=[0.01,0.1,0.5,0.7,0.9,0.99]
e_L=[0.01,0.1,0.2,0.3,0.5,0.7,0.8,0.9,0.99]
#e_L=[0.01,0.1]
#a_L=[0.01,0.1,,0.3,0.5,0.8,1.2,0.5,]#a=0.001
XI=[]
for k in range(1,d):
XI.append(f.fftfreq(2*N)*2*N*2*np.pi/(2*L))
XI.append(f.rfftfreq(2*N)*(2*N)*2*np.pi/(2*L))
xi=np.array(np.meshgrid(*XI,indexing='ij'))
#k2=xi**2+eta**2
X= np.mgrid[tuple(slice(- L, L, dx) for _ in range(d))].T
Y=X.reshape((-1,d))
I=np.ones_like(Y.sum(axis=-1))
if corr=="norm2":
p=2
elif corr=="norm1":
p=1
AX=[-k for k in range(1,d+1)]
dist=sp.distance_matrix(Y,Y,p=p)
samples=6 #several disorder per value of epsilon to reduce the error
nb_eq_moy=[]
nb_eq_std=[]
nb_eq_tot=[]
points_meta=[]
for e in e_L:
v=1-e
a=e
def gamma_V_func(r,corr=corr):
if corr=="norm2":
return v**2*np.exp(-1/2*(r/l_cV)**2)
elif corr=="norm1":
return v**2*np.exp(-r/l_cV)
def gamma_A_func(r,corr=corr):
if corr=="norm2":
return a**2*np.exp(-1/2*(r/l_cA)**2)
elif corr=="norm1":
return a**2*np.exp(-r/l_cA)
eps=np.identity(len(Y))*10**-12
def make_V(size=1,d=d,ampl=v,l_c=l_cV,xi=xi,N=N,dist=dist,eps=eps):
covV=gamma_V_func(dist)
L=cholesky(covV+eps,lower=True)
if corr=="norm2":
k2=(xi**2).sum(axis=0)
covVH=(2*N)**d*ampl*(2*np.pi*l_c**2)**(d/2)*np.exp(-l_c**2*k2/4)
elif corr=="norm1":
truc=2*l_c/(1+(l_c*xi)**2)
covVH=(2*N)**d*ampl*2*truc.prod(axis=0)
if size>1:
listV=[]
listKV=[]
for k in range(size):
Vh=covVH*np.random.standard_normal(covVH.shape)
V=1/2*f.irfftn(Vh)
V2=V.reshape(Y[:,0].shape)
KV=cho_solve((L,True),V2)
listV.append(V)
listKV.append(KV)
return listV,listKV
else:
Vh=covVH*np.random.standard_normal(covVH.shape)
V=1/2*f.irfftn(Vh)
V2=V.reshape(Y[:,0].shape)
KV=cho_solve((L,True),V2)
return V,KV
def make_A(size=1,d=d,ampl=a,l_c=l_cA,xi=xi,N=N,dist=dist,eps=eps,corr=corr):
covA=gamma_A_func(dist)
L=cholesky(covA+eps,lower=True)
A=np.zeros((d,d)+(2*N,2*N))
KA=np.zeros((d,d,(2*N)**d))
if corr=="norm2":
k2=(xi**2).sum(axis=0)
covAH=(2*N)**d*ampl*(2*np.pi*l_c**2)**(d/2)*np.exp(-l_c**2*k2/4)
elif corr=="norm1":
truc=2*l_c/(1+(l_c*xi)**2)
covAH=(2*N)**d*ampl*2*truc.prod(axis=0)
if size>1:
list_A=[]
list_KA=[]
for k in range(size):
for i in range(1,d):
for j in range(i):
AH=covAH*np.random.standard_normal(covAH.shape)
A[i,j,:]=1/2*f.irfftn(AH)
A[j,i,:]=-A[i,j,:]
A2=A[i,j,:].reshape(Y[:,0].shape)
KA[i,j,:]=cho_solve((L,True),A2)
KA[j,i,:]=-KA[i,j,:]
list_A.append(A)
list_KA.append(KA)
return list_A,list_KA
else:
for i in range(1,d):
for j in range(i):
AH=covAH*np.random.standard_normal(covAH.shape)
A[i,j,:]=1/2*f.irfftn(AH)
A[j,i,:]=-A[i,j,:]
A2=A[i,j,:].reshape(Y[:,0].shape)
KA[i,j,:]=cho_solve((L,True),A2)
KA[j,i,:]=-KA[i,j,:]
return A,KA
def Gamma_VX(x,Y,corr=corr):
x=x.reshape((-1,1,x.shape[-1]))
if corr=="norm2":
r=npa.sqrt(((Y-x)**2).sum(axis=-1))
return v**2*npa.exp(-1/2*(r/l_cV)**2)
elif corr=="norm1":
r=npa.abs(Y-x).sum(axis=-1)
return v**2*npa.exp(-r/l_cV)
def Gamma_AX(x,Y):
x=x.reshape((-1,1,x.shape[-1]))
if corr=="norm2":
r=npa.sqrt(((Y-x)**2).sum(axis=-1))
return a**2*npa.exp(-1/2*(r/l_cA)**2)
elif corr=="norm1":
r=npa.abs(Y-x).sum(axis=-1)
return a**2*npa.exp(-r/l_cA)
Vl,KVl=make_V(samples)
Al,KAl=make_A(samples)
points_L=[]
nb_eq_L=[]
nb_L=[]
for p in range(samples):
V=Vl[p]
A=Al[p]
KV=KVl[p]
KA=KAl[p]
def V_func(x,Y=Y):
return npa.dot(Gamma_VX(x,Y),KV)
def A_func(x,Y=Y):
return npa.tensordot(Gamma_AX(x,Y),KA,[[1],[2]])
nablaV=egrad(V_func,argnum=0)
coordA={}
for i in range(1,d):
for j in range(i):
def A_coord(x,Y=Y):
AT=npa.tensordot(Gamma_AX(x,Y),KA,[[1],[2]]).T
return AT[i,j,:]
coordA[(i,j)]=egrad(A_coord)
def jacobA(x,coordA=coordA):
JA=np.zeros((d,d)+x.shape)
for i in range(1,d):
for j in range(i):
JA[i,j,:]=coordA[(i,j)](x)
JA[j,i,:]=-JA[i,j,:]
return JA
def rotA(x):
jacA=jacobA(x)*1/np.sqrt(2)
return np.einsum('ij...j -> ...i',jacA)
Time=7*10**2
dt=2*10**-3
T=0.0
threshold=dt*(T+0.1)
def run(j=1,Time=Time,dt=dt,T=T,mu=mu,threshold=threshold):
np.random.seed()
x=np.random.uniform(-1,1,size=2)
points=np.empty((Time,2))
fixe=0
K=0
for k in range(Time):
points[k]=x
u=-nablaV(x)+np.random.standard_normal(size=2)*T+rotA(x)#-mu*x
x_new=x+dt*u
x[0]=(1+x_new[0])%2-1
x[1]=(1+x_new[1])%2-1
if k>50 and np.std(points[k-20:k,0])+np.std(points[k-20:k,1])<threshold: #stop condition: stop trajectories which have converged to something
fixe=1
K=k
return [x[0],x[1],fixe,K]
Nt=200
th=24
if __name__=='__main__':
out=np.array(process_map(run,range(Nt),max_workers=th)).T #this one is parallelized on realizations
else:
points=dist(Nt)
x=out[0]
y=out[1]
fixed=out[2]
K=out[3]
points=np.array([x,y]).T
points_eq=points[fixed==1]
points_eq=points_eq.T
points_L.append(points)
hi=np.histogram2d(points_eq[0],points_eq[1],bins=20)
nb=hi[0]
nb_L.append(nb)
nb_eq=len(nb[nb>0])
nb_eq_L.append(nb_eq)
nb_L=np.array(nb_L)
nb_eq_L=np.array(nb_eq_L)
points=np.array(points_L)
points_meta.append(points_L)
nb_eq_tot.append(nb_eq_L)
nb_eq_moy.append(nb_eq_L.mean())
nb_eq_std.append(nb_eq_L.std())
params=np.array([d,L,N,l_cV,l_cA,v,a,mu,Time,dt])
#data=[V,KV,A,KA,points,params]
#np.savez('points_per_norm1_2.npz',e_L=e_L,nb_eq_std=nb_eq_std,nb_eq_moy=nb_eq_moy,points=points_meta,params=params)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment