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

start

parents
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
d=2
#dx=0.02
L=1
N=32
dx=L/N
#N=L/dx
mu=0
l_cV=0.3
l_cA=0.3
alpha_A=1
alpha_V=1
v=1.0
a=1.0
def gamma_V_func(r):
return v**2*np.exp(-1/2*(r/l_cV)**(2*alpha_V))
def gamma_A_func(r):
return a**2*np.exp(-1/2*(r/l_cA)**(2*alpha_A))
xi,eta=f.rfftfreq(2*N)*(2*N)*2*np.pi/(2*L), f.fftfreq(2*N)*2*N*2*np.pi/(2*L)
xi, eta=np.meshgrid(xi,eta)
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))
dist=sp.distance_matrix(Y,Y,p=2)
covA=gamma_A_func(dist)
eps=np.identity(len(Y))*10**-12
#f.rfft2(np.random.standard_normal(X[:,:,0].shape))
def make_V(d=d,ampl=v,l_c=l_cV,k2=k2,N=N,dist=dist,eps=eps):
covV=gamma_V_func(dist)
L=cholesky(covV+eps,lower=True)
Vh=(2*N)**d*ampl*(2*np.pi*l_c**2)**(d/2)*np.exp(-l_c**2*k2/4)*np.random.standard_normal(k2.shape)
V=1/2*f.irfft2(Vh)
V2=V.reshape(Y[:,0].shape)
KV=cho_solve((L,True),V2)
return V,KV
def make_A(d=d,ampl=a,l_c=l_cA,k2=k2,N=N,dist=dist,eps=eps):
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))
for i in range(1,d):
for j in range(i):
AH=(2*N)**d*ampl*(2*np.pi*l_c**2)**(d/2)*np.exp(-l_c**2*k2/4)*np.random.standard_normal(k2.shape)
A[i,j,:]=1/2*f.irfft2(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()
def Gamma_VX(x,Y):
x=x.reshape((-1,1,x.shape[-1]))
r=npa.sqrt(((Y-x)**2).sum(axis=-1))
return v**2*npa.exp(-1/2*(r/l_cV)**(2*alpha_V))
def Gamma_AX(x,Y):
x=x.reshape((-1,1,x.shape[-1]))
r=npa.sqrt(((Y-x)**2).sum(axis=-1))
return a**2*npa.exp(-(r/l_cA)**alpha_A)
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]])
internship_random_force @ 9c37b47a
Subproject commit 9c37b47afb08cc1ef82083a51320e01794ce3fcc
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, det
from numpy import fft as f
from tqdm import tqdm
import multiprocessing as mp
from tqdm.contrib.concurrent import process_map
from jax.ops import index, index_add, index_update
from jax.config import config
import scipy.optimize as opti
import time
import finufft
import os
import warnings
#warnings.simplefilter("ignore",np.ComplexWarning)
os.environ['KMP_DUPLICATE_LIB_OK']='True'
config.update('jax_platform_name', 'cpu')
config.update("jax_enable_x64", True)
d=2
L=np.pi
N=64
dx=L/N
mu=0
corr_L=["norm1","norm2"]
corr="norm2"
l_cV=0.14*np.pi
l_cA=0.14*np.pi
#e_L=[0.01,0.5,0.9,0.99]
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(d):
XI.append(f.fftfreq(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
AX=[-k for k in range(1,d+1)]
samples=6
th=24
#Z=[Y.T[k] for k in range(d)]
#Z=tuple(Z)
def make_V(size=1,d=d,l_c=l_cV,xi=xi,N=N,corr=corr,nb=len(e_L)):
if corr=="norm2":
k2=(xi**2).sum(axis=0)
covVH=(2*np.pi*l_c**2)**(d/2)*np.exp(-l_c**2*k2/4) #(2*N)**d*
elif corr=="norm1":
truc=2*l_c/(1+(l_c*xi)**2)
covVH=2*truc.prod(axis=0) #(2*N)**d*
if size>1:
list_Vh=[]
for k in range(size):
Vh=covVH*(np.random.standard_normal(covVH.shape)+1j*np.random.standard_normal(covVH.shape))
list_Vh.append(Vh)
list_Vh=np.array(list_Vh)
return list_Vh
else:
Vh=covVH*(np.random.standard_normal(covVH.shape)+1j*np.random.standard_normal(covVH.shape))
return Vh
def make_A(size=1,d=d,l_c=l_cA,xi=xi,N=N,corr=corr,nb=len(e_L)):
Ah=np.zeros((d,d)+xi[0].shape,dtype=np.complex128)
if corr=="norm2":
k2=(xi**2).sum(axis=0)
covAH=(2*np.pi*l_c**2)**(d/2)*np.exp(-l_c**2*k2/4) #(2*N)**d*
elif corr=="norm1":
truc=2*l_c/(1+(l_c*xi)**2)
covAH=2*truc.prod(axis=0) #(2*N)**d
if size>1:
list_Ah=[]
for k in range(size):
for i in range(1,d):
for j in range(i):
Ah[i,j,:]=covAH*(np.random.standard_normal(covAH.shape)+1j*np.random.standard_normal(covAH.shape))
Ah[j,i,:]=-Ah[i,j,:]
list_Ah.append(Ah)
list_Ah=np.array(list_Ah)
return list_Ah
else:
for i in range(1,d):
for j in range(i):
Ah[i,j,:]=covAH*(np.random.standard_normal(covAH.shape)+1j*np.random.standard_normal(covAH.shape))
Ah[j,i,:]=-Ah[i,j,:]
return Ah
Vh=make_V()
Ah=make_A()
e=0.5
v=1-e
a=e
nabVh=v*np.array([1j*xi[k]*Vh for k in range(d)])
rotAh=a*np.array([np.array([1j*xi[j]*Ah[i,j,:] for j in range(d)]).sum(axis=0) for i in range(d)])/np.sqrt(d)
##
N2=20
dx2=L/N2
X2= np.mgrid[tuple(slice(- L, L, dx2) for _ in range(d))].T
Y2=X2.reshape((-1,d))
#Zt=[Y2.T[k] for k in range(d)]
#Z=tuple(Zt)
fh=-nabVh+rotAh
fhx=np.array([1j*xi[k]*fh[0] for k in range(d)])
fhy=np.array([1j*xi[k]*fh[1] for k in range(d)])
plan=finufft.Plan(nufft_type=2,n_modes_or_dim=(2*N,)*d,modeord=1,eps=1e-10,isign=1,n_trans=d,fftw=64)
#plan.setpts(*Z)
#t1=time.time()
#force=np.real(plan.execute(fh).reshape(X2.shape))
def force2(x,fh=fh,plan=plan):
plan.setpts(*(x.reshape((d,1))))
force=np.real(plan.execute(fh))
return (force**2).sum(axis=-1)
def Dforce2(x,fh=fh,fhx=fhx,fhy=fhy,plan=plan):
plan.setpts(*(x.reshape((d,1))))
force=np.real(plan.execute(fh))
fx=np.real(plan.execute(fhx))
fy=np.real(plan.execute(fhy))
return force[0]*fx+force[1]*fy
#print((force2<=10**-4).sum()/4)
#ppoints=X2[(force2<=10**-4)].reshape((-1,d))
#hi=np.histogramdd(ppoints,bins=200)
#nb=hi[0]
#nb_eq=len(nb[nb>0])
minim=[]
for y in Y2:
res=opti.minimize(force2,y,jac=Dforce2,bounds=[(-L,L)]*d)
x=res.x
plan.setpts(*(x.reshape((d,1))))
force=np.real(plan.execute(fh))
if np.sqrt((force**2).sum())<=10**-5:
minim.append(x)
minim=np.array(minim).reshape((-1,d))
stable=[]
for x in minim:
jac=np.zeros((2,2))
plan.setpts(*(x.reshape((d,1))))
jac[0]=plan.execute(fhx)
jac[1]=plan.execute(fhy)
if det(jac)>0 and np.trace(jac)<0:
stable.append(x)
stable=np.array(stable).reshape((-1,d))
hi=np.histogramdd(stable,bins=100)
nb=hi[0]
nb_eq=len(nb[nb>0])
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