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

rework

parent 16e3a7e8
File added
# internship_random_force
Repository containing all relevant codes resulting from my internship at Econophysix during spring 2021 on the equilibria of random forces field in low dimension.
\ No newline at end of file
import numpy as np
import jax.numpy as jnp
from jax import grad, jit, vmap,jacfwd
from jax import random
import scipy.spatial as sp
from scipy.linalg import cholesky, cho_solve
from numpy import fft as f
from jax.config import config
config.update('jax_platform_name', 'cpu')
config.update("jax_enable_x64", True)
d=2
L=4
N=64
dx=L/N
mu=0
corr_L=["norm1","norm2","long"]
#corr="norm2"
corr="long"
gamma=0.8
t1=time.time()
l_cV=0.05
l_cA=0.05
e=0.0
#e_L=[0.01,0.1,0.2,0.3,0.5,0.7,0.8,0.9,0.99]
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'))
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
elif corr=="long":
p=2
AX=[-k for k in range(1,d+1)]
dist=sp.distance_matrix(Y,Y,p=p)
samples=6
th=24
eps=np.identity(len(Y))*10**-12
def gamma_V_func(r,corr=corr,gamma=gamma):
if corr=="norm2":
return np.exp(-1/2*(r/l_cV)**2)
elif corr=="norm1":
return np.exp(-r/l_cV)
elif corr=="long":
return (1+l_cV**-2*r**2)**(-gamma/2)
def gamma_A_func(r,corr=corr,gamma=gamma):
if corr=="norm2":
return np.exp(-1/2*(r/l_cA)**2)
elif corr=="norm1":
return np.exp(-r/l_cA)
elif corr=="long":
return (1+l_cA**-2*r**2)**(-gamma/2)
def Gamma_VX(x,Y,corr=corr,gamma=gamma):
if corr=="norm2":
r=jnp.sqrt(((Y-x)**2).sum(axis=-1))
return jnp.exp(-1/2*(r/l_cV)**2)
elif corr=="norm1":
r=jnp.abs(Y-x).sum(axis=-1)
return jnp.exp(-r/l_cV)
elif corr=="long":
r=jnp.sqrt(((Y-x)**2).sum(axis=-1))
return (1+l_cV**-2*r**2)**(-gamma/2)
def Gamma_AX(x,Y):
if corr=="norm2":
r=jnp.sqrt(((Y-x)**2).sum(axis=-1))
return jnp.exp(-1/2*(r/l_cA)**2)
elif corr=="norm1":
r=jnp.abs(Y-x).sum(axis=-1)
return jnp.exp(-r/l_cA)
elif corr=="long":
r=jnp.sqrt(((Y-x)**2).sum(axis=-1))
return (1+l_cA**-2*r**2)**(-gamma/2)
def make_V(size=1,d=d,l_c=l_cV,xi=xi,N=N,corr=corr,eps=eps,gamma=gamma):
L=cholesky(gamma_V_func(dist)+eps,lower=True)
if corr=="norm2":
k2=(xi**2).sum(axis=0)
covVH=(2*N)**d*(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*2*truc.prod(axis=0)
elif corr=="long":
k=l_c*np.sqrt((xi**2).sum(axis=0))+10**-9
bd=-(gamma-d)/2
covVH=(2*N)**d*2*np.pi**(d/2)*(k/2)**bd*kv(bd,k)/Gamma(bd+1)
if size>1:
list_V=[]
list_KV=[]
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)
list_V.append(V)
list_KV.append(KV)
list_V=np.array(list_V)
list_KV=np.array(list_KV)
return list_V,list_KV
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,l_c=l_cA,xi=xi,N=N,eps=eps,corr=corr,gamma=gamma):
L=cholesky(gamma_A_func(dist)+eps,lower=True)
A=np.zeros((d,d)+(2*N,)*d)
KA=np.zeros((d,d,(2*N)**d))
if corr=="norm2":
k2=(xi**2).sum(axis=0)
covAH=(2*N)**d*(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*2*truc.prod(axis=0)
elif corr=="long":
k=np.sqrt((xi**2).sum(axis=0))+10**-9
bd=(d-gamma)/2
covAH=(2*N)**d*2*np.pi**(d/2)*(k/2)**bd*kv(bd,k)/Gamma(bd+1)
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)
listA=np.array(list_A)
listKA=np.array(list_KA)
return listA,listKA
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
V,KV=make_V()
A,KA=make_A()
v=1-e
a=e
def V_func(x,Y=Y,v=v):
return v*jnp.dot(Gamma_VX(x,Y),KV)
nablaV=vmap(grad(V_func,argnums=0))
def A_func(x,Y=Y,a=a):
return a*(Gamma_AX(x,Y)*KA).sum(axis=-1)
jacobA=jacfwd(A_func)
cste=1/np.sqrt(d)
def rotA(x,d=d,cste=cste):
jacA=jacobA(x)*cste
return jnp.einsum('ij...j -> ...i',jacA)
rotA=vmap(rotA)
V_f=vmap(V_func)
V2=V_f(Y).reshape(V.shape)
print((V2-V).std())
\ No newline at end of file
internship_random_force @ 9c37b47a
Subproject commit 9c37b47afb08cc1ef82083a51320e01794ce3fcc
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