In [6]:
from pylab import *
%matplotlib inline
for iq in range(shape(Pval)[0]):
plot(Omega, Pval[iq,:].real/nF)
plot(Omega, Pval[iq,:].imag/nF)
show()
First we copy things that we developed in the class.
import time
from scipy import *
from mweight import *
from numpy import linalg
from numpy import random
from scipy import special
import sys
import numpy as np
from numba import jit
from pylab import *
from numpy import random
class params:
def __init__(self):
self.kF = 1. # typical momentum
self.cutoff = 3*self.kF # integration cutoff
self.dkF = 0.1*self.kF # the size of a step
#
self.Nitt = 2000000 # total number of MC steps
self.Ncout = 50000 # how often to print
self.Nwarm = 1000 # warmup steps
self.tmeassure = 10 # how often to meassure
#
self.Nbin = 129 # how many bins for saving the histogram
self.V0norm = 2e-2 # starting V0
self.dexp = 6 # parameter for fm at the first iteration, we will use
self.recomputew = 2e4/self.tmeassure # how often to check if V0 is correct
self.per_recompute = 7 # how often to recompute fm auxiliary measuring function
@jit(nopython=True)
def TrialStep0(qx, momentum):
tiQ = int( random.rand()*len(qx) ) # trial iQ for qx[iQ]
Ka_new = qx[tiQ] # |r_0| length of the vector
th, phi = pi*random.rand(), 2*pi*random.rand() # spherical angles for vector q in spherical coordinates
sin_th = sin(th) # trial step probability is proportional to sin(theta) when using sperical coodinates
Q_sin_th = Ka_new * sin_th
K_new = array([Q_sin_th*cos(phi), Q_sin_th*sin(phi), Ka_new*cos(th)]) # new 3D vector r_0
q2_sin2_old = sum(momentum[0,:2]**2) # x^2+y^2 = r_old^2 sin^2(theta_old)
q2_old = q2_sin2_old + momentum[0,2]**2 # x^2+y^2+z^2 = r_old^2
trial_ratio = 1.
if q2_old != 0: # make sure we do not get nan
sin_th_old = sqrt(q2_sin2_old/q2_old) # sin(theta_old)
if sin_th_old != 0: # make sure we do not get nan
trial_ratio = sin_th/sin_th_old
accept=True
return (K_new, Ka_new, trial_ratio, accept, tiQ)
@jit(nopython=True)
def TrialStep1(iloop,momentum,dkF,cutoff):
dk = (2*random.rand(3)-1)*dkF # change of k in cartesian coordinates of size p.dkF
K_new = momentum[iloop,:] + dk # K_new = K_old + dK
Ka_new = linalg.norm(K_new) # norm of the new vector
trial_ratio = 1. # trial step probability is just unity
accept = Ka_new <= cutoff # we might reject the step if the point is too far from the origin
return (K_new, Ka_new, trial_ratio, accept)
@jit(nopython=True)
def Give_new_X(momentum, K_new, iloop):
tmomentum = copy(momentum)
tmomentum[iloop,:] = K_new # this is trial configuration X_{new}=momentum
return tmomentum
This is implementation of the Lindhardt formula:
\begin{eqnarray} P(q,\omega) = -2 \int \frac{d^3k\; d\Omega_q}{(2\pi)^3} \frac{f(\varepsilon_{\vec{k}+\vec{q}})-f(\varepsilon_{\vec{k}})}{\omega-\varepsilon_{\vec{k}+\vec{q}}+\varepsilon_\vec{k}+i\delta} \end{eqnarray}where $\varepsilon_\vec{k}=k^2-k_F^2$.
The input is momentum, where momentum[0,:] is $\vec{q}$ and momentum[1,:] is $\vec{k}$. It returns the value for array of frequencies.
@jit(nopython=True)
def ferm(x):
if x>700:
return 0.
else:
return 1./(exp(x)+1.)
@jit(nopython=True)
def Linhard2_inside(momentum, Omega, kF, T, broad, nrm):
q = momentum[0,:] # Notice that q is the first component of momenta
k = momentum[1,:]
e_k_q = linalg.norm(k-q)**2 - kF*kF
e_k = linalg.norm(k)**2 - kF*kF
dfermi = (ferm(e_k_q/T)-ferm(e_k/T))
return -2*nrm*dfermi/(Omega-e_k_q+e_k+broad*1j)
class Linhard2:
def __init__(self, Omega, kF, T, broad):
self.Omega = Omega
self.kF = kF
self.T = T
self.broad = broad
self.nrm = 1/(2*pi)**3
self.Ndim = 2 # Notice we need to add this for metropolis routine
def __call__(self, momentum):
return Linhard2_inside(momentum, self.Omega, self.kF, self.T, self.broad, self.nrm)
Now we modify IntegrateByMetropolis so that func is allowed to return complex array of values.
We marked all lines that require change by
# Next line needs CORRECTION for homework
Pval[iq,iOm] becomes $q$ and $\Omega$ array.
When we compute ratio, we need to take only the zero frequency value of the function, i.e., abs(fQ_new[0][0]) and fQ[0][0].
In measurement, the weight W should depend only zero frequency value only, to cancel the probability from the Metropolis, i.e., abs(fQ[0][0])
When measuring, we need to properly update the values Pval[iQ,:]+= f0.
The physical weight $\widetilde{V}_{physical}$ for determining if V0 is of the right magnitude, needs to be changed so that it depends on zero frequency only f0[0]
When we reset the stored values in Pval, we need to create array of proper size and type
Pval = zeros(shape(Pval), dtype=complex)
Printing should be modified, for example ratio has to be correctly modified to depend on zero frequency obly, and only the zero frequency value of the function is printed out.
def IntegrateByMetropolis4(func, qx, p):
""" Integration by Metropolis:
func(momentum) -- function to integrate
qx -- mesh given by a user
p -- other parameters
Output:
Pval(qx)
"""
tm1 = time.time()
random.seed(0) # make sure that we always get the same sequence of steps
Pval = zeros((len(qx),len(func.Omega)),dtype=complex) # CHANGE: Final results V_physical is stored in Pval
Pnorm = 0.0 # V_alternative is stored in Pnorm
Pval_sum = 0.0 # this is widetilde{V_physical}
Pnorm_sum = 0.0 # this is widetilde{V_alternative}
V0norm = p.V0norm # this is V0
dk_hist = 1.0 # we are creating histogram by adding each configuration with weight 1.
Ndim = func.Ndim # dimensions of the problem
inc_recompute = (p.per_recompute+0.52)/p.per_recompute # How often to self-consistently recompute
# the wight functions g_i and h_{ij}.
momentum = zeros((Ndim,3)) # contains all variables (r1,r2,r3,....r_Ndim)
# We call them momentum here, but could be real space vectors or momentum space vectors.
iQ = int(len(qx)*random.rand()) # which bin do we currently visit for r0, iQ is current r0=qx[iQ]
momentum[1:,:] = random.random((Ndim-1,3)) * p.kF / sqrt(3.) # Initial guess for r1,r2,....r_N is random
momentum[0,:] = [0,0,qx[iQ]] # initial configuration for r_0 has to be consistent with iQ, and will be in z-direction
# This is fm function, which is defined in mweight.py module
mweight = meassureWeight(p.dexp, p.cutoff, p.kF, p.Nbin, Ndim) # measuring function fm in alternative space
# fQ on the current configuration. Has two components (f(X), V0*f_m(X))
fQ = func(momentum), V0norm * mweight( momentum ) # fQ=(f(X), V0*f_m(X))
#print('starting with f=', fQ, '\nstarting momenta=', momentum)
t_sim, t_mes, t_prn, t_rec = 0.,0.,0.,0.
Nmeassure = 0 # How many measurements we had?
Nall_q, Nall_k, Nall_w, Nacc_q, Nacc_k = 0, 0, 0, 0, 0
c_recompute = 0 # when to recompute the auxiliary function?
for itt in range(p.Nitt): # long loop
t0 = time.time()
iloop = int( Ndim * random.rand() ) # which variable to change, iloop=0 changes external r_0
accept = False
if (iloop == 0): # changing external variable : r_0==Q
Nall_q += 1 # how many steps changig external variable
(K_new, Ka_new, trial_ratio, accept, tiQ) = TrialStep0(qx, momentum)
else: # changing momentum ik>0
Nall_k += 1 # how many steps of this type
(K_new, Ka_new, trial_ratio, accept) = TrialStep1(iloop,momentum,p.dkF,p.cutoff)
if (accept): # trial step successful. We did not yet accept, just the trial step.
tmomentum = Give_new_X(momentum, K_new, iloop)
fQ_new = func(tmomentum), V0norm * mweight(tmomentum) # f_new
# Notice that we take |f_new(X)+V0*fm_new(X)|/|f_old(X)+V0*fm_old(X)| * trial_ratio
ratio = (abs(fQ_new[0][0])+fQ_new[1])/(abs(fQ[0][0])+fQ[1]) * trial_ratio # !!CHANGE!!
accept = abs(ratio) > 1-random.rand() # Metropolis
if accept: # the step succeeded
momentum[iloop] = K_new
fQ = fQ_new
if iloop==0:
Nacc_q += 1 # how many accepted steps of this type
iQ = tiQ # the new external variable index
else:
Nacc_k += 1 # how many accepted steps of this type
t1 = time.time()
t_sim += t1-t0
if (itt >= p.Nwarm and itt % p.tmeassure==0): # below is measuring every p.tmeassure steps
Nmeassure += 1 # new meassurements
W = abs(fQ[0][0])+fQ[1] # !!CHANGE!! this is the weight we are using
f0, f1 = fQ[0]/W, fQ[1]/W # the two measuring quantities
Pval[iQ,:]+= f0 # CHANGE: V_physical : integral up to a constant
Pnorm += f1 # V_alternative : the normalization for the integral
Pnorm_sum += f1 # widetilde{V}_alternative, accumulated over all steps
Wphs = abs(f0[0]) # !!CHANGE!! widetilde{V}_{physical}, accumulated over all steps
Pval_sum += Wphs
# doing histogram of the simulation in terms of V_physical only.
# While the probability for a configuration is proportional to f(X)+V0*fm(X), the histogram for
# constructing g_i and h_{ij} is obtained from f(X) only.
mweight.Add_to_K_histogram(dk_hist*Wphs, momentum, p.cutoff, p.cutoff)
if itt>10000 and itt % (p.recomputew*p.tmeassure) == 0 :
# Now we want to check if we should recompute g_i and h_{ij}
# P_v_P is V_physical/V_alternative*0.1
P_v_P = Pval_sum/Pnorm_sum * 0.1
# We expect V_physical/V_alternative*0.1=P_v_P to be of the order of 1.
# We do not want to change V0 too much, only if P_V_P falls utside the
# range [0.25,4], we should correct V0.
change_V0 = 0
if P_v_P < 0.25 and itt < 0.2*p.Nitt: # But P_v_P above 0.25 is fine
change_V0 = -1 # V0 should be reduced
V0norm /= 2 # V0 is reduced by factor 2
Pnorm /= 2 # V_alternative is proportional to V0, hence needs to be reduced too.
Pnorm_sum /= 2 # widetilde{V}_alternative also needs to be reduced
if P_v_P > 4.0 and itt < 0.2*p.Nitt: # and P_v_P below 4 is also fine
change_V0 = 1 # V0 should be increased
V0norm *= 2 # actually increasing V0
Pnorm *= 2
Pnorm_sum *= 2
if change_V0: # V0 was changed. Report that.
schange = ["V0 reduced to ", "V0 increased to"]
print('%9.2fM P_v_P=%10.6f' % (itt/1e6, P_v_P), schange[int( (change_V0+1)/2 )], V0norm )
# Here we decied to drop all prior measurements if V0 is changed.
# We could keep them, but the convergence can be better when we drop them.
Pval = zeros(shape(Pval), dtype=complex) # !!CHANGE!!
Pnorm = 0
Nmeasure = 0
# Next we should check if g_i and h_ij need to be recomputed.
# This should not be done too often, and only in the first half of the sampling.
if (c_recompute==0 and itt<0.7*p.Nitt):
t5 = time.time()
# At the beginning we recompute quite often, later not so often anymore
# as the per_recompute is increasing...
p.per_recompute = int(p.per_recompute*inc_recompute+0.5)
# We normalized f_m, hence all previous accumulated values are now of the order
# of 1/norm. We also normalize the new additions to histogram with similar value,
# but 5-times larger than before.
dk_hist *= 5*mweight.Normalize_K_histogram()
if dk_hist < 1e-8: # Once dk becomes too small, just start accumulating with weight 1.
dk_hist = 1.0
mweight.Recompute()# Here we actually recompute g_i and h_{ij}.
fQ = func(momentum), V0norm * mweight( momentum ) # And now we must recompute V0*f_m, because f_m has changed!
#print(shape(fQ[0]),fQ[1])
t6 = time.time()
print('%9.2fM recomputing f_m=%10.6f' % (itt/1e6, fQ[1]))
t_rec += t6-t5
c_recompute += 1
if c_recompute>=p.per_recompute : c_recompute = 0 # counting when we will recompute next.
t2 = time.time()
t_mes += t2-t1
if (itt+1)% p.Ncout == 0 : # This is just printing information
P_v_P = Pval_sum/Pnorm_sum * 0.1 # what is curent P_v_P
Qa = qx[iQ] # current r0
ka = linalg.norm(momentum[1,:]) # current r1
ratio = (abs(fQ_new[0][0])+fQ_new[1])/(abs(fQ[0][0])+fQ[1]) # CHANGE. current ratio
print( '%9.2fM Q=%5.3f k=%5.3f fQ_new=%8.3g,%8.3g fQ_old=%8.3g,%8.3g P_v_P=%10.6f' % (itt/1e6, Qa, ka, abs(fQ_new[0][0]), fQ_new[1], abs(fQ[0][0]), fQ[1], P_v_P) ) # CHANGE
t3 = time.time()
t_prn += t3-t2
Pval *= len(qx) * V0norm / Pnorm # Finally, the integral is I = V0 *V_physical/V_alternative
# This would be true if we are returning one single value. But we are sampling len(qx) values
# And we jump between qx[i] uniformly, hence each value should be normalized with len(qx).
tp1 = time.time()
print('Total acceptance rate=', (Nacc_k+Nacc_q)/(p.Nitt+0.0), 'k-acceptance=', Nacc_k/(Nall_k+0.0), 'q-acceptance=', Nacc_q/(Nall_q+0.0))
print('k-trials=', Nall_k/(p.Nitt+0.0), 'q-trial=', Nall_q/(p.Nitt+0.0) )
print('t_simulate=%6.2f t_meassure=%6.2f t_recompute=%6.2f t_print=%6.2f t_total=%6.2f' % (t_sim, t_mes, t_rec, t_prn, tp1-tm1))
return (Pval,mweight)
from numpy import random
rs = 2.
kF = pow( 9*pi/4., 1./3.) / rs # given in the homework, this is kF of the uniform electron gas
nF = kF/(2*pi*pi) # this is density of states at the fermi level in uniform electron gas
T = 0.02*kF**2 # temperature, as given in the homework
broad = 0.002*kF**2 # broadening, as given in the homework
cutoff = 3*kF # cutoff for the integration, as given in the homework
Omega = linspace(0,kF**2,100)
qx = linspace( 0.1*kF, 0.4*kF, 4)
lh2 = Linhard2(Omega, kF, T, broad)
p = params()
p.Nitt = 5000000
(Pval,mweight) = IntegrateByMetropolis4(lh2, qx, p)
0.02M recomputing f_m= 0.005468
0.04M P_v_P= 0.248214 V0 reduced to 0.01
0.05M Q=0.096 k=0.966 fQ_new= 0.0326, 0.00432 fQ_old= 0.0352, 0.00601 P_v_P= 0.494056
0.10M Q=0.384 k=0.940 fQ_new= 0.0888, 0.00596 fQ_old= 0.0888, 0.00596 P_v_P= 0.474734
0.15M Q=0.384 k=0.944 fQ_new= 0.0393, 0.00596 fQ_old= 0.0393, 0.00596 P_v_P= 0.474189
0.18M recomputing f_m= 0.004361
0.20M Q=0.192 k=0.984 fQ_new= 0.0804, 0.00436 fQ_old= 0.0804, 0.00436 P_v_P= 0.476963
0.25M Q=0.288 k=0.990 fQ_new=5.41e-06,0.000528 fQ_old= 0.0717, 0.00436 P_v_P= 0.475829
0.30M Q=0.192 k=0.977 fQ_new= 0.0264, 0.00548 fQ_old= 0.0264, 0.00548 P_v_P= 0.474379
0.35M Q=0.384 k=0.943 fQ_new= 0.0778, 0.00531 fQ_old= 0.0778, 0.00531 P_v_P= 0.473089
0.36M recomputing f_m= 0.000782
0.40M Q=0.288 k=1.058 fQ_new= 1.6e-06, 0.00112 fQ_old= 1.6e-06, 0.00112 P_v_P= 0.472408
0.45M Q=0.096 k=1.015 fQ_new=0.000557, 0.00279 fQ_old= 0.0139, 0.00279 P_v_P= 0.474629
0.50M Q=0.384 k=1.037 fQ_new= 0.0342, 0.00173 fQ_old= 0.0342, 0.00173 P_v_P= 0.474981
0.55M Q=0.096 k=0.931 fQ_new= 0.0124, 0.00538 fQ_old= 0.0694, 0.00538 P_v_P= 0.475035
0.56M recomputing f_m= 0.005537
0.60M Q=0.288 k=0.949 fQ_new= 0.00168, 0.0027 fQ_old= 0.102, 0.00543 P_v_P= 0.476021
0.65M Q=0.096 k=0.994 fQ_new= 0.00117, 0.00439 fQ_old= 0.00117, 0.00439 P_v_P= 0.477640
0.70M Q=0.192 k=0.931 fQ_new= 0.0199, 0.00543 fQ_old= 0.0199, 0.00543 P_v_P= 0.478103
0.75M Q=0.096 k=0.975 fQ_new= 0.00631, 0.00554 fQ_old= 0.0505, 0.00554 P_v_P= 0.477862
0.78M recomputing f_m= 0.001137
0.80M Q=0.096 k=0.958 fQ_new= 0.103, 0.0055 fQ_old= 0.103, 0.0055 P_v_P= 0.478747
0.85M Q=0.192 k=0.961 fQ_new= 0.0346, 0.0055 fQ_old= 0.104, 0.0055 P_v_P= 0.478020
0.90M Q=0.096 k=0.978 fQ_new= 0.0559, 0.00438 fQ_old= 0.0556, 0.00438 P_v_P= 0.478010
0.95M Q=0.288 k=0.797 fQ_new= 0.0192,0.000804 fQ_old= 0.0192,0.000804 P_v_P= 0.476874
1.00M Q=0.384 k=0.830 fQ_new= 0.0267, 0.00105 fQ_old= 0.0267, 0.00105 P_v_P= 0.476792
1.02M recomputing f_m= 0.000106
1.05M Q=0.384 k=0.654 fQ_new= 0.0138,0.000164 fQ_old= 0.0138,0.000164 P_v_P= 0.476931
1.10M Q=0.096 k=1.274 fQ_new= 3.8e-20,3.13e-05 fQ_old=1.62e-18,4.87e-05 P_v_P= 0.477011
1.15M Q=0.384 k=1.062 fQ_new= 0.0128, 0.00115 fQ_old= 0.0128, 0.00115 P_v_P= 0.477801
1.20M Q=0.096 k=0.946 fQ_new= 0.0536, 0.00532 fQ_old= 0.0536, 0.00532 P_v_P= 0.478066
1.25M Q=0.192 k=0.889 fQ_new=7.19e-05, 0.00276 fQ_old= 0.0364, 0.00276 P_v_P= 0.479229
1.28M recomputing f_m= 0.002761
1.30M Q=0.096 k=0.969 fQ_new= 0.0617, 0.00432 fQ_old= 0.0987, 0.00543 P_v_P= 0.478654
1.35M Q=0.192 k=1.037 fQ_new= 0.0246, 0.0017 fQ_old= 0.0246, 0.0017 P_v_P= 0.479435
1.40M Q=0.288 k=0.951 fQ_new=7.02e-06, 0.00113 fQ_old= 0.104, 0.00533 P_v_P= 0.479198
1.45M Q=0.384 k=1.052 fQ_new=6.15e-07, 0.00113 fQ_old=6.15e-07, 0.00113 P_v_P= 0.479911
1.50M Q=0.192 k=1.024 fQ_new= 0.0278, 0.0017 fQ_old= 0.0278, 0.0017 P_v_P= 0.478871
1.55M Q=0.192 k=1.018 fQ_new=1.79e-05, 0.00276 fQ_old= 0.0375, 0.00276 P_v_P= 0.478718
1.56M recomputing f_m= 0.005446
1.60M Q=0.288 k=0.998 fQ_new=0.000745, 0.0041 fQ_old= 0.0264, 0.00437 P_v_P= 0.478599
1.65M Q=0.288 k=0.960 fQ_new= 0.00998, 0.00545 fQ_old= 0.00998, 0.00545 P_v_P= 0.478511
1.70M Q=0.096 k=0.970 fQ_new= 0.0221, 0.00545 fQ_old= 0.102, 0.00545 P_v_P= 0.478637
1.75M Q=0.288 k=0.723 fQ_new=1.81e-10,0.000361 fQ_old= 0.0166,0.000361 P_v_P= 0.478983
1.80M Q=0.192 k=0.901 fQ_new=9.48e-05, 0.00278 fQ_old=9.48e-05, 0.00278 P_v_P= 0.479206
1.85M Q=0.096 k=0.902 fQ_new= 0.0423, 0.00278 fQ_old= 0.0423, 0.00278 P_v_P= 0.479127
1.86M recomputing f_m= 0.005472
1.90M Q=0.384 k=0.834 fQ_new= 0.026, 0.00105 fQ_old= 0.026, 0.00105 P_v_P= 0.479072
1.95M Q=0.288 k=1.008 fQ_new= 0.0268, 0.00275 fQ_old= 0.0268, 0.00275 P_v_P= 0.478795
2.00M Q=0.288 k=1.036 fQ_new=2.66e-05, 0.00173 fQ_old= 0.0213, 0.00173 P_v_P= 0.478583
2.05M Q=0.288 k=1.117 fQ_new= 6.3e-10,0.000383 fQ_old= 0.0151,0.000383 P_v_P= 0.478373
2.10M Q=0.096 k=0.955 fQ_new= 0.0359, 0.00547 fQ_old= 0.0877, 0.00547 P_v_P= 0.478332
2.15M Q=0.192 k=0.849 fQ_new= 0.00994, 0.00144 fQ_old= 0.0233, 0.00144 P_v_P= 0.478204
2.18M recomputing f_m= 0.001392
2.20M Q=0.384 k=0.940 fQ_new= 0.0113, 0.00536 fQ_old= 0.0113, 0.00536 P_v_P= 0.478434
2.25M Q=0.192 k=0.986 fQ_new= 0.0018, 0.00437 fQ_old= 0.07, 0.00437 P_v_P= 0.478629
2.30M Q=0.192 k=0.954 fQ_new= 0.0471, 0.00544 fQ_old= 0.0471, 0.00544 P_v_P= 0.478481
2.35M Q=0.096 k=0.996 fQ_new= 0.00534, 0.00437 fQ_old= 0.00534, 0.00437 P_v_P= 0.478497
2.40M Q=0.288 k=0.899 fQ_new= 0.0167, 0.0028 fQ_old= 0.0167, 0.0028 P_v_P= 0.478139
2.45M Q=0.192 k=1.110 fQ_new= 0.00309,0.000411 fQ_old= 0.0157,0.000576 P_v_P= 0.478104
2.50M Q=0.096 k=0.885 fQ_new= 0.0253, 0.0028 fQ_old= 0.0331, 0.0028 P_v_P= 0.478316
2.52M recomputing f_m= 0.004406
2.55M Q=0.288 k=0.847 fQ_new= 0.0282, 0.0014 fQ_old= 0.00106, 0.0014 P_v_P= 0.478575
2.60M Q=0.096 k=0.955 fQ_new= 0.0259, 0.00548 fQ_old= 0.0279, 0.00548 P_v_P= 0.478553
2.65M Q=0.288 k=0.966 fQ_new= 0.0312, 0.00548 fQ_old= 0.0312, 0.00548 P_v_P= 0.478621
2.70M Q=0.288 k=1.072 fQ_new=7.43e-07,0.000793 fQ_old= 0.0205,0.000793 P_v_P= 0.478200
2.75M Q=0.384 k=0.838 fQ_new= 0.0119, 0.0014 fQ_old= 0.0119, 0.0014 P_v_P= 0.478403
2.80M Q=0.192 k=1.096 fQ_new= 0.0127,0.000577 fQ_old= 0.0127,0.000577 P_v_P= 0.478360
2.85M Q=0.192 k=1.014 fQ_new=6.98e-05, 0.00277 fQ_old=6.98e-05, 0.00277 P_v_P= 0.478273
2.88M recomputing f_m= 0.004189
2.90M Q=0.096 k=0.921 fQ_new= 0.00133, 0.00419 fQ_old= 0.00133, 0.00419 P_v_P= 0.478150
2.95M Q=0.096 k=0.972 fQ_new= 0.00678, 0.00559 fQ_old= 0.0319, 0.00559 P_v_P= 0.478048
3.00M Q=0.288 k=0.747 fQ_new=1.24e-08,0.000472 fQ_old= 0.018,0.000472 P_v_P= 0.478427
3.05M Q=0.192 k=0.976 fQ_new= 0.087, 0.00559 fQ_old= 0.087, 0.00559 P_v_P= 0.478643
3.10M Q=0.096 k=0.931 fQ_new= 0.0585, 0.00549 fQ_old= 0.0585, 0.00549 P_v_P= 0.479055
3.15M Q=0.384 k=0.954 fQ_new=5.97e-05, 0.0028 fQ_old= 0.0211, 0.00559 P_v_P= 0.479024
3.20M Q=0.096 k=0.661 fQ_new=8.64e-13,0.000134 fQ_old=8.64e-13,0.000134 P_v_P= 0.478824
3.25M Q=0.192 k=0.993 fQ_new=5.67e-06,0.000744 fQ_old= 0.0355, 0.00446 P_v_P= 0.479098
3.26M recomputing f_m= 0.004406
3.30M Q=0.096 k=0.971 fQ_new= 0.0231, 0.00549 fQ_old= 0.0231, 0.00549 P_v_P= 0.479112
3.35M Q=0.384 k=1.026 fQ_new= 0.0174, 0.00169 fQ_old= 0.0275, 0.00169 P_v_P= 0.479120
3.40M Q=0.288 k=0.897 fQ_new= 0.00175, 0.00285 fQ_old= 0.0165, 0.00285 P_v_P= 0.479208
3.45M Q=0.288 k=0.944 fQ_new= 0.0223, 0.00541 fQ_old= 0.0223, 0.00541 P_v_P= 0.478889
3.50M Q=0.384 k=1.046 fQ_new= 0.016, 0.00169 fQ_old= 0.016, 0.00169 P_v_P= 0.478991
3.55M Q=0.096 k=0.981 fQ_new= 0.0746, 0.00441 fQ_old= 0.0746, 0.00441 P_v_P= 0.478867
3.60M Q=0.384 k=1.181 fQ_new=1.49e-11,0.000193 fQ_old= 0.0128,0.000193 P_v_P= 0.478942
3.65M Q=0.384 k=0.888 fQ_new= 0.00513, 0.00108 fQ_old= 0.0333, 0.00285 P_v_P= 0.478775
3.70M Q=0.384 k=0.951 fQ_new= 0.0145, 0.00541 fQ_old= 0.0145, 0.00541 P_v_P= 0.478789
3.75M Q=0.192 k=1.136 fQ_new=6.18e-05,0.000381 fQ_old=3.22e-11,0.000381 P_v_P= 0.478741
3.80M Q=0.096 k=0.921 fQ_new=0.000156,0.000843 fQ_old= 0.0446, 0.00416 P_v_P= 0.479042
3.85M Q=0.192 k=0.934 fQ_new= 0.0388, 0.00541 fQ_old= 0.0388, 0.00541 P_v_P= 0.479213
3.90M Q=0.288 k=0.918 fQ_new=0.000491, 0.00416 fQ_old=0.000491, 0.00416 P_v_P= 0.479223
3.95M Q=0.384 k=0.985 fQ_new= 0.0306, 0.00441 fQ_old= 0.0306, 0.00441 P_v_P= 0.479165
4.00M Q=0.288 k=0.966 fQ_new=9.31e-09,0.000544 fQ_old= 0.0201, 0.00549 P_v_P= 0.478880
4.05M Q=0.288 k=1.073 fQ_new= 0.0184,0.000753 fQ_old= 0.0184,0.000753 P_v_P= 0.479174
4.10M Q=0.288 k=1.145 fQ_new=9.66e-12,0.000273 fQ_old= 0.0172,0.000273 P_v_P= 0.479025
4.15M Q=0.192 k=0.818 fQ_new=2.06e-05, 0.00108 fQ_old=2.06e-05, 0.00108 P_v_P= 0.478914
4.20M Q=0.192 k=0.942 fQ_new= 0.0688, 0.00549 fQ_old= 0.0694, 0.00541 P_v_P= 0.478924
4.25M Q=0.288 k=0.963 fQ_new= 0.0181, 0.00549 fQ_old= 0.0927, 0.00549 P_v_P= 0.478838
4.30M Q=0.288 k=1.034 fQ_new= 0.0193, 0.00169 fQ_old= 0.0193, 0.00169 P_v_P= 0.478744
4.35M Q=0.288 k=0.853 fQ_new= 0.015, 0.00142 fQ_old= 0.015, 0.00142 P_v_P= 0.478821
4.40M Q=0.096 k=0.932 fQ_new= 0.0014, 0.00541 fQ_old= 0.0525, 0.00541 P_v_P= 0.478944
4.45M Q=0.288 k=1.012 fQ_new= 0.0261, 0.00277 fQ_old= 0.0261, 0.00277 P_v_P= 0.478897
4.50M Q=0.096 k=0.992 fQ_new=7.09e-08,0.000753 fQ_old= 0.00776, 0.00441 P_v_P= 0.479101
4.55M Q=0.288 k=0.778 fQ_new=1.28e-08,0.000664 fQ_old= 0.0182,0.000664 P_v_P= 0.479017
4.60M Q=0.096 k=0.944 fQ_new= 0.0144, 0.00541 fQ_old= 0.0906, 0.00541 P_v_P= 0.478976
4.65M Q=0.096 k=0.937 fQ_new= 0.00175, 0.00142 fQ_old= 0.0436, 0.00541 P_v_P= 0.478913
4.70M Q=0.384 k=0.932 fQ_new= 0.0234, 0.00541 fQ_old= 0.02, 0.00541 P_v_P= 0.478956
4.75M Q=0.096 k=0.849 fQ_new=8.55e-06, 0.00142 fQ_old=8.55e-06, 0.00142 P_v_P= 0.479104
4.80M Q=0.192 k=0.922 fQ_new=0.000712, 0.00416 fQ_old= 0.0602, 0.00416 P_v_P= 0.479204
4.85M Q=0.288 k=0.891 fQ_new= 0.0168, 0.00285 fQ_old= 0.0168, 0.00285 P_v_P= 0.479194
4.90M Q=0.096 k=0.987 fQ_new= 0.00405, 0.00441 fQ_old= 0.00405, 0.00441 P_v_P= 0.479026
4.95M Q=0.192 k=0.966 fQ_new= 0.0281, 0.00549 fQ_old= 0.0281, 0.00549 P_v_P= 0.478891
5.00M Q=0.192 k=1.049 fQ_new= 0.0205, 0.00109 fQ_old= 0.0205, 0.00109 P_v_P= 0.479008
Total acceptance rate= 0.4979084 k-acceptance= 0.5911432016364302 q-acceptance= 0.4046314467069558
k-trials= 0.500113 q-trial= 0.499887
t_simulate= 54.79 t_meassure= 20.17 t_recompute= 0.02 t_print= 1.31 t_total= 78.10
from pylab import *
%matplotlib inline
for iq in range(shape(Pval)[0]):
plot(Omega, Pval[iq,:].real/nF)
plot(Omega, Pval[iq,:].imag/nF)
show()